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ficient algorithms for implementing collocation on the whole sphere are described. 

A formal relationship between collocation and least squares adjustment is used to ob¬ 
tain an alternative form of the collocation algorithm that is likely to be stable with 
dense data sets and, with a minor modification, can be used to implement least 
squares adjustment as well. The basic principle is that for regular grids the variance* 
covariance matrix of the data consists of Toeplitz-circulant blocks, so it can be both 
set up and inverted very efficiently. 

Section four describes a method for computing the covariances between area 
means that is much more efficient than the usual approach of integrating numerically 
the point covariance function twice over the blocks. This method Is essential for 
setting up the variance-covariance matrix of die data to implement collocation when the 
data are area means. In this and other respects the present work is complementary 
to Report 291, by the same author. 
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Finally, there is an Appendix containing the listing of subroutines that imple¬ 
ment various algorithms, together with some explanatory notes on how to use them. 
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1. Introduction 


Spherical harmonics are closely associated with the basic theory of 
gravitational and magnetic fields, such as those of the Earth and planets; 
for this reason they are important both in geodesy and in Earth and planetary 
physics. 

The present work considers the numerical aspects of the reduction of global 
data sets to spherical harmonic coefficients, so the emphasis has been laid on the 
algorithms for this purpose. The procedures for harmonic analysis (and synthesis! 
given here are general enough to be used in the study of magnetic or electric fields, 
but most conclusions regarding their accuracy are restricted to the gravity field 
of our planet and to fields with the same power spectrum. The accuracy of these 
methods cannot be separated from the type of signal being used. 

Modern instrumentation has provided scientists and engineers with vast 
amounts of information, and modern computers have made the processing of it 
possible, and even routine, thanks to constant improvements in both hardware 
and software. In the mid sixties, those branches of applied mathematics, physics, 
and engineering concerned with the sifting of data, or with the study of very large 
regular structures, were greatly affected by the advent of the Fast Four.er Trans¬ 
form (FFT). In spite of the fact that spherical harmonics are members of the 
family of Fourier transforms, closely related to two dimensional Fourier series, 
geodesy has lagged behind in the development of techniques similar to the FFT, 
partly because of the rather wicked nature of the sphere on which data are usually 
given, partly because there have not been enough data to make the development of 
powerful techniques a general concern. The topological differences between the 
Euclidean plane and the surface of the sphere may very well prevent the finding 
of algorithms as efficient as the FFT for the latter (certainly none seems to have 
been reported to date) but such algorithms should be regarded, nonetheless, as 
the desideratum for all those who wish to put their time and work into developing 
good numerical methods for spherical harmonic analysis. 

The increasing use of artificial satellites for surveying the gravity field, 
particularly by radar altimetry and by the projected tracking of one satellite by 
another, are making the use of very efficient and accurate techniques for handling 
the resulting data imprescindible; even a casual review of the literature of the 
last few years will show that serious efforts to provide such techniques are getting 
under way. The days of scarce, scattered, unreliable data are just about over. 

The remainder of this section defines the basic problem and the associated 
notation, shows the relationships between spherical harmonics and 2D-Fourier 
series, presents some of the similitudes and differences between both, and ex¬ 
plains some efficient algorithms for harmonic analysis and synthesis that are 
common to a number of different problems. Section 2 begins by defining a 
quadratic measure for the accuracy of the estimated harmonic coefficients based 
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on the covariance functions of the signal and the noise; a discussion on the 
optimization of this measure follows, leading to the application of least squares 
collocation to harmonic analysis. The use and implementation of least squares 
adjustment follows, and then a discussion of the connection between least squares 
and least squares collocation , shown as alternative and efficient techniques 
for solving the same problem. The section closes with an algorithm for the case 
when data are Irregularly distributed. Section 3 illustrates with several 
numerical examples some of the methods presented earlier. Section 4 
introduces an efficient formula for computing the covariances between block 
averages when analyzing mean values by collocation. This formula Is much 
more efficient than others based on the numerical quadratures of the "point" 
covariance function. 


1.1. Spherical Harmonic Analysis and Synthesis; Definitions 


A square integrable, analytical function f(0,A) defined on the unit sphere 
0 £ 0 £ tt and 0 £ X < 2 tt can be expanded in a series of surface spherical 
harmonics 


f(0,A) 



L Pn» (cos 0) [ C RB cos mX + S ni 


■ = 0 


sinmA ] 


( 1 . 1 ) 


where: P nI are the associated Legendre functions of the first kind, fully 

normalized so-i—. f P B , (cos 0) 3 f cos> \ mA da = 1 (here f da 
4 rr Jo \Sin/ Jo 

indicates integration on the unit sphere); 

C at ,S n , are the fully normalized spherical harmonic coefficients. 


For the sake of brevity, the following alternative notation shall be used when possible 


-a a , _ rP„(cos 0)cos mA if a = 0 

, a » ’ lp ni (cos 0) sinmA if a = 1 

and 

« * _ f C* if Oi a 0 

u ‘ " is„. if a=i 

~“ /y 

The purpose of spherical harmonic analysis is to estimate the coefficients C a , 
from measurements of the signal f ( 0, A). These measurements, which may be 
corrupted by some noise or error signal "n", and which are assumed in what 
follows to be finite in number, constitute the data. The individual samples are 
called z t j, so z tJ = f(0 t , A,) + n tJ . The subscripts i and j are used only to 
designate the position of the sample in a two-dimensional array, or grid, covering 
in some more or less regular way the sphere: i corresponds always to (co)latitude, 
and j to longitude. While, as in the last paragraph of section 3, some places in 
the grid may be empty, the grid itself is defined by a set of complete parallels 
and meridians . Unless otherwise specified, the separation between the lines of 
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latitude can be variable, but that between meridians is always constant and 
equal to AX = t/N, where N is an integer. The grid most often consi¬ 
dered in this work is the equal angular grid, also called the regular grid; 
in this case &B is also constant and equal to A A . 

For equal angular grids i and j take values 0 £ i £ N-l and 0 ^ j £2N-1, 
i increasing from North to South, and j from West to East. Formulas where the 
i,j subscripts appear are in the form appropriate to the regular grid, though most 
of them can be extended in a very simple way to other partitions. 

Data may consist of values determined at the intersections of the grid, in 
which case they are referred to as "point data ", or they may be averages over 
the blocks defined by the lines of the grid, and then they are called "area means " 
or 1 'block means ". In the equal angular grid the northernmost and southernmost 
blocks reach to the respective poles; there are N rows of blocks (i.e. blocks be¬ 
tween the same parallels), there are 2N blocks per row, and i,j identify 
olocks according to the row and column they are in. Grids for equal angular 
point data may be "center point" grids, data being measured at the center of 
each block, so i = 0 corresponds to 9 = A = AX/2. 


No blocks stride the equator, so regular grids are symmetrical with respect 
to it; in other words; N is always even (extension to N odd is trivial). 


Area means are identified here by overbars; they can be expanded in series 
simply by integrating (1.1) term by term, which can be done because the spherical 
harmonic series is always uniformly convergent for 0 s 0 5 n and 0 £ X £ 2 -r ; 

fu = “ E J. J Y“(9,A)da 

i 1 


'! 3 n = 0 « = 0 O =° 

os r. 1 


= —- E E E P a9 Pr.«(cosS)sinSd8 [c n . cosmA 
Ajj „ =0 .so er=o J e t L J X. 

+ S BB ^sinmXdxJ (1.2) 

Here Ti 3 is the area mean of f( 0,X) on the block a, 3 whose area is 
A[j = AX(cos - cos (9. + A0)) • 


1 if at = 3, m = p and n = k 
otherwise 


A basic property of spherical harmonics is their orthogonality: 

{J 

as a consequence of which 

C® = jirJo ?*<0,A) *<3,A) da 


(1.3) 


(1.4) 
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Expression (1.4) is the inverse of (1.1); both constitute! the basis of spherical 
harmonic analysis. In general, (1.4) cannot be calculated analytically, because 
what is available is not the function f (0,A), but a finite set of noisy measure¬ 
ments in the form of point values z u or their area means zj 3 . Discretizing 
(1.4) on an equal angular grid results, for instance, in the following numerical 
quadratures formula; 


A 

C 


a 

n» 


1 

4tt 


N-l 3N-1 

Z Z Y«(St.*j) f( Oi,X 4 ) 

1 = 0 j = o 


A.j 


(1.5) 


a a _a 

where C n « indicates the estimate of C na . as this type of formula is usually only 
an approximation. Formulas resembling (1.5) can be called "point valuer-type 
quadrature formulas',’ and can be handled by algorithms that are all identical from 
a structural point of view and whose prototype is that of paragraph 1.5. 


If the data are block averages, several simple approximations have been 
proposed that take the form 

A a N-l dN-l _ 

C«> - /i, E Z ^ij J Y na (0, A) do 
1=0 1=0 

where the are scale factors. This kind of formula shall be studied further 
in section 2. Writing the above expression in full 

r 8,*^8 f V«rco S1 

C„. = Mn ^ 0 f u J fll P„.(COS0) sin9d0J^ isin/ mXdA < 1 * G ) 
This belongs to the general type 

N-1 SN -1 

& “ 1 E r,I [U!m!} C0Smi ' lX + {A!mS} SiDmiW ] (1 - 7 > 


where B(m) = | 


f(cosm AA-l)/m if m^O 


if m =0 ’ A < m > = l 


_ f(sin m AA/ m if m / 0 


AA 


if m = 0 


Formulas resembling (1.6) or (1. 7) appear several times in this work, and are 
called here area means-type quadrature formulas. 


-.O' 

If all C 


the (N max + if terms in 

f( 0 ,X) w max 


»« with 0 ^ n ^ N max 


^ Z E Z c 

n = 0 " = 0 Oi — O 


are known, they can be used to compute 
“ Y.® ( 0, A) d.8) 


which can be regarded as an approximation to f ( 0 , A) at the point ( 9 , A) . 
Expression (1.8), together with the truncation at N max of (1.2) (area means), 
defines the object of spherical harmonic synthesis ; given the coefficients, 
estimate the function. As shown later, analysis and synthesis are related by 
a simple duality, so they require about the same number of operations when 
performed on a certain grid and with a given number of coefficients. 
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(1.9) 


The set of all degree variances 

= i c„, 2 + s a / 

• = o 

constitutes the power spectrum of f(0,X). If all coefficients of degree n > N m ax 
are zero, so are the corresponding <J n ; insuchcase the function f( 0,X) is said 
to be band limited. 


Closely associated with the power spectrum is the isotropic covariance 
function 

CO 

cov(f(P),f(Q) ) = E a n s P n (cos ((>«,) (1.10) 

0 

Here p and Q are two points on the sphere, tfa = cos^f cos 9? cos 0 Q sin 6 p 
x sinSij cos( Ap - Xq) ] is the geocentric angle or spherical distance between both, 
and P B (cos<^) is the unnormalized Legendre polynomial of degree n. The power 
spectrum and the covariance of a function on the sphere , like those of func¬ 
tions defined on the real line, on the plane,and on Euclidean spaces of higher 
dimensions, can be used to obtain estimators related to discrete Wiener 
filters and predictors, i. e., satisfying a minimum variance condition. The 
spherical version of such technique is known as least squares collocation (V.ontz, 
1972). The inverse of expression (1.10) is 

a a 3 = (2n + l) J o P„(cos lilpq) cov(f(P),f(Q) ) sin d^ (1.11) 

Equations (1.10) and (1.11) show that, as usual, power spectrum and covariance 
function are linear transforms of each other. A formula such as (1.11) is some¬ 
times called a Legendre transform. 


In addition to the autocovariance (1.10), one can define, more generally, 

CO 

COV(U(P),V(Q) ) = Yj C n ,V Pn(COS tfpq) ( 1 . 10 *) 

n ss 0 

as the covariance of two functions u(0p, Xp) and v( d Q , X q ) on the sphere, where 
the n 


c«"’ = E Cn a , c; +S B U . s n v . , n = 0,1,... (1.9 

«=o 


U,v 


*\ 


constitute the "power crosspectrum”. 

A relationship similar to (1.11) applies to the c Y and to cov (u(P), v(Q)). 

A very important property of spherical harmonics is Parseval's theorem: 

-tluB.Xf da - E ( 1 . 12 , 
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The orthogonality of spherical harmonics, and the fact that they form a 
complete orthogonal set of functions on the sphere, are among the reasons why 
they are used so widely in theory and in practice, but there is more to them than 
orthogonality. 


The solid spherical harmonics 


-prn-?“(0.A) and r a Y®(9,X) 


where r is the distance to the origin of coordinates, are all solutions of Laplace's 
equation,which in Cartesian coordinates is 


v a w 


_ a 3 w ^ a 3 w . 

“ *r + *?- 


5 3 W 

3z a 


= 0 


(1. 13) 


This makes them appropriate for the study of harmonic functions, such as the 
gravitational potential, in spherical coordinates. Another property they have, 
unique among functions that are orthogonal on the sphere, is the relationship 


[* f (P) E (2n + l)k* P,(0«) da = E E E k »C?.Y!(Q) (1-14) 

4 TT JC7 l = 0 n -r-Q q ss 0 QL =0 

which corresponds to the Convolution Theorem for ordinary Fourier series* and is 
the basis of such fundamental formulas as Stokes’ in gravimetric geodesy. 


1.2. Relationship Between Spherical Harmonics and 2-D Fourier Series 


As indicated in the previous paragraph, spherical harmonics share important 
properties with ordinary (trigonometric) Fourier series in one or more dimensions. 
There is also a very immediate relationship with ordinary two-dimensional (2-D) 
series that will be explained here. From Hobson (1931, Ch. HI, formula (7)) we 
know that 


P a „(cos Q) 


- (-D" (2n)t 

2* n! (n- m)! 


Sin-8 (cos—8- ttjag .Ija .- J) 
l 2 (2n -1) 


cos 


9 + 


so the normalized 


(n- m)., ■ (n - m - 3) 
+ 2* 4 • (2n-l)(2n- 3) 

Legendre function is 


cos 


a- » “4 i 


(1.15) 


P B1 (cos 9) 

where 

L(n, m) = 


(-1) (2n!) 
2* n! (n-m)! 

(n-m)/2 if 
(n-m-l)/2 if 


^/tSEESjEEIsto-e t ) a.(n.m ) cos—“a 


n-m is even 
n-m is odd 


while a„ (n, m) 


(-l)* (n-m)(n-m-l). ..(n-m-2k + l)[2 • 4...2k(2n-1)...(2n-2k+l)] _l 

for k > 0 
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In the interval -n £ 9 ^ 71 sin* 9 is even when m is even, and odd when 
m is odd. In the interval -ns 9 £ n the sum L ^^a k (n m) cos “~ 8 '" 8 ''9 

k = o 

is always even (it is a sum of even powers of cos 9 ), so the parity of 
P a (cos 9) Is the same as that of m if -tt s e s -. An even function can 
be expanded into a sum of cosines. and an odd function into a sum of sines . 
The highest frequency term will correspond to the highest frequency in the 
expansion of sin® 9 cos’*® 9, so this term will be of the form a* cos n9 or 
b» sinn9. Therefore, the Legendre function satisfies one of the following 
equations in -tt £ 9 £ n: 

a 

Pea (cos 9 ) = £ C\ cos t0 (a) if m is even ; 

t =0 

__ a 

Pna (cos 9) = £ S* t * sint9 (b) if m is odd 

t = o 

-a 

A spherical harmonic Y» (9,A) can have one of four possible forms: 


for m even < 


—Of 


<9, A) * 


for m odd: 

_a 

Y«. 


n 

z 

, t=o 


(9 A) - 


Xo 

a 

z 

(t = 0 


C. cos tQ 

cos mX 

for 

a = 0 

C* B eos te 

sinmX 

for 

II 

0 

C** sint0 

cos mA 

for 

a = 0 

n n 

C t sint9 

sinmX 

for 

p 

11 

M 


(1.16a) 


(1.16b) 


A sum of spherical harmonics such as (1. 8 ) is equivalent to a sum of terms 
of the form cf t B sint 9 cosmA, C* B sintg sinmX, C* n cost9 cosmA and 
(f t n costs sinmX, which are also the basic functions of 2-D Fourier series. 

The highest m and n in the spherical harmonic expansion 

"•« a I _ _ a 

f (9,X) = [ [ [ Ca. Y«. (9,X) 

are equal to Na. » . so the highest degree and order (or spatial frequencies) 
in the Fourier series are also equal to Nn, x . In conclusion: every surface 
spherical harmonic expansion where the highest degree is N® ax is identical 
to a 2-D Fourier series (where the highest m and n are also Nb„), 
in the domain -ti M s n , 0 s X s 2n . The converse is not true, be¬ 
cause continuous functions on a sphere, such as the Y bD , must satisfy certain 
conditions at the poles that ordinary functions on the -rrsgsrr, 0 s X s 2 n 
domain do not have to. Spherical harmonics correspond to a subclass 
(linear subspace) of 2-D Fourier series. 






For example: 


(Heiskanen and Moritz, Chapter 1,19G7) 
P n (cose) = /l sin0 

(cose) =/30 sing cos8 = -§^0 sin26 


So 

Y a ° = /3 sine cosX, Y u x =/3 sinG sinX 
Ya!= 1^30 sin 20 cos2X, Yz* = %/3 0 sin20 sin2X 
Calling the 2 -D Fourier coefficients "a P „",where 


a£. 

correspond 

to 

terms 

of the 

form 

cospG cos raX 


it 

tt 

If 

t! I! 

ft 

cospG sinmX 

a?« 

ir 

»i 

ft 

If If 

II 

sinpG cos mX 

4. 

it 

it 

tt 

tt tt 

tf 

sinp9 sinmX 


these can be related to the respective Cb» by expressions of the form 

C®, « T JU^\)&= M± £ r* >p af« n=m, m + 1, m + 

» (n+m) ! p-o 


2 , . . . 

(1.17) 


where (3 = a if m is even, and j3 = 0£+ 2 if m is odd. The l“ f p are 
defined as 


T a 


n, p 


In 


. P 



cos p0 P n , (cos 8) sin G d 9 
sinpG ft B (cos 6) sin G d0 


if m is even 


if m is odd. 


and can be computed recursively using the formula 
2n - 


U 


*"- p ~ 2(n 
with the following starting values 


- 1 . V* ) n + m - 1 

- m) I • p+ j n - m 


Tn- a, p 


( o if (m + p) is odd 


L a, p 


(1.18a) 


(1.18b) 


2(m+ 1) (2m) ! 


if m is 


2" [(m + l) 3 - p 3 ] [(m - l) 3 - p 3 ]. . .[3 s - p 2 ] [l a - p a ] even, p even, 

_ 2 fm-t- 1^ 12m) i __ if m is odd, 

2*0 m+ l) 3 - 1?] [(m - l) 2 - p 2 ] . . . \f - p 3 ] [-p] p odd. 


The l" &p are zero for alternate values of both p and m. These equations 
were reported by Ricardi and Burrows (1972), and show how to obtain the 





— CX ft 

C« a once the a have been computed from the data by means of the 2-D 
discrete Fourier transform. Normally this would be impossible, because 
the data only exists in the upper half of the -it s 9 i , 0 £ X ^ 2rr 
interval, and the_3-D algorithm requires information on all of it. However, 
the fact that the P„, (cos 9) are sums of sines only or cosines onl}' makes 
the calculation possible. 

While the maximum degree and order in a 2-D Fourier series d o not 
reach the Nyquist frequency 1 (m, n < N = in an equal angular grid, 
all the coefficients can be recovered exactly by solving (2N) 2 equations 
such as 


UN 3 

EE z 

lao i=9 B = o 


JCOS 

jsin 



mX, 


(1.19) 


When n or m exceed N, iX , the matrix of the system of e quations 
becomes singular, and the discrete Fourier trans form consists 
of coefficients that "fit" the data, but differ from the true coefficients. The 
estimated coefficients are said to have been aliased with those that exceed 
the Nyquist frequency. In the case of spherical harmonics , which are a 
special case of 2-D Fourier series, a similar situation must arise: the har¬ 
monics in the data with n s N are going to be aliased with those of lower 
degree, so the information available is not enough to recover all coefficients 
because the sampling is too coarse . 

The aliasing of spherical harmonics sampled on regular grids is a 
consequence of the aliasing of the respective 2-D Fourier series, so it 
makes sense to talk of a "Nyquist frequency" in the case of those functions. 
Having established the connection between aliasing in both types of series, 
it is time to point out also some Important differences. 


1.3 Sampling Errors 

Expressions (1.16a b) shows that spherical harmonics are fini te 
sums of 2-D Fourier harmonics, which is not the same as being each a 
Fourier harmonic. From this simple fact follow some important dis¬ 
tinctions. 


To understand them better, let us begin by stating some basic pro¬ 
perties of Fourier series in one dimension, which carry over to higher 
dimensions but are easier to explain in one dimension. 

~T 

If sampled at a constant interval AX = — , the following is always true 
of sines and of cosines : * 


1 Named after the Nyquist Theorem: 
of period 2NAX can be recovered only if 


the Fourier coefficients of a function 
N**x < N. 
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if m ^ p < N 


3N-1 


cos mkAX cos pkA X = 0 


a n-\ 

£ sinmkAX sinpkAX = 0 
v = o 
9lU 

]T cos mkAX sinpkAX = 0 

k = 0 

Ttese expressions are discrete counterparts of 

r 3TT /- an 

/ cosmX cospXdX = / sinmX sinpXdX = 0 

Jo J 0 ...1 __ 


if m ^p < N 

for all m and all p. 


(]. 20-a) 


jTT 


f Q cosmX sinpXdX=0 


when m ,^p 
for all m and p. 


(1.20-b) 


and show that the orthogonal properties of sines and cosines are maintained 
when these are sampled regularly, provided the Nyquist frequency is not 
exceeded. From this follows that 

3 N » 1 3 TT 

a^ = H T (sin) rnk AXf(kAX) = hf (sin) mXf( >) dX (1.21) 

y=o J 0 


where H = 


4 if m = 0 

I otherwise ^ h = 


in if m = 0 

- otherwise 
rr 


, so the "Fourier” 


counterpart of (1.5), (i.e. (1.21)) is an exact "numerical quadratures" 
formula for Fourier series. When the Nyquist frequency is exceeded, the 
trigonometric relationships 


cos mkAX= cos (2N t m) kA X 
- sinmkaX=sin(2N * m)kAX 


imply that (1.21) will give, not the true coefficients, but the aliased estimates 

K 


a°, = a° 


+ E a anN+« + E 

h=0 h=1 

fll « = a *» + E al ahN+» E 


a 0 

a 3hN-a 


( 1.2 2 ) 


h =o 


h =al 


2 hN-a 


with K such that 2KN is below, and 2(K+1)N is above the highest frequency present in 
f(X). Expression (1. 20 a-b), (1. 21), and (1.22) are the foundations of discrete 
Fourier analysis (also known as the computation of the discrete Fourier 
transform, or D. F.T.), and so well known that they are almost second nature 
to many engineers and scientists. Unfortunately, none of these discrete 
formulas has exact counterparts in spherical harmonic analysis, and this 
fact has been the cause of considerable confusion. The most common 
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misunderstandings seem to occur in the following areas: 


(a) Number of equations and number of unknowns : 

At each point in the grid it is possible to write an equation such as 

N Ba x n I Ct ft 

f(C'iX 1 )= £ £ £ c„ Y n , (B u Xi) (1.23) 

a = o » = o Q| = 0 

similar to (1.19). The maximum number of coefficients that can be rcco- 
vered vvitliout exceeding die Nyquist frequency (i. e., those with n < N) 
is N 2 . Coefficients of degree or order equal or higher than N shall not 

a , 

be, in general, free from aliasing. .Since there are 2N~ points in an 
equiangular grid, it follows that the maximum cumber of fully recoverable 
coefficients is also half the number of points (equations) in the grid. By con¬ 
trast, in 2HD Fourier analysis the number of coefficients equals the number 
of points in the grid (in the interval si < 6 s n, 0 < 1 -< 2 ti). One could 
have a grid with only N 2 points, as proposed by Giacaglia and Lundquist 
(1972), but such a grid would not be equal angular on the sphere. 


The author had discussed this problem elsewhere (Colombo, 

1979a, paragraph (5.2) ), showing that the system of equations (1.23) 
becomes singular when all coefficients of degree and order 0 e (» , m) s m 
are included among the unknowns, and M £ N. In other words: it is not 
possible to solve for a complete set of coefficients to degree and order M ^ N. 
The relevant part of that argument can be summarized as follows: the 
columns of A f the_giatrix of the system of equations (1.23), consist in succes¬ 
sions of values Y n „ ( e lf A.,) of the harmonics corresponding t o the un¬ 
knowns C®„ at the points (0^ Aj) in the grid. The scalar product of two 
such columns is. 


N-12N-1 _ N _ _ 

£ £ y“ Y* = £ P BB (cose 5 ) Pp< (COS0,) 

1= o i _ 0 lr:0 


3N-1 

* z 


3 = ° 


( COS ) . ) cos) 

) sin (rnjAX ) sin = 


0 if o 1? p 
0 If m / q 
7*0 otherwise 


aocording to (1.20-a). Therefore, if two columns correspond to unknowns 
of different orders m and q, they must be orthogonal and, thus, indepen¬ 
dent. For the whole matrix not to be singular, all columns of the same order 
m must form sub-matrices A(m) that have full rank. Otherwise, there 
will be columns in those A(m), and consequently in A, that are linearly 
dependent, so A cannot bo inverted. Consider A(0), corresponding to all 
unknowns of order 0. This is a 2N s x (M + l)jmatrix, and the elements of 
the columns of A(0) have the same values as P n0 (cos0,), 0 -s n £ M, 

0 ^ i * N-l. The P n „ are functions of e* only, and there are N par¬ 
allels in the grid, so there are nc. more than N independent rows in A(0). 
Because the P r „(cosO) are polynomials of degree n in oos G, these 
rows form a sub-matrix S(0) of A(0) that has M + 1 independent columns. 







as long asM + l^N or M < N. If M a N the extra columns will turn S(0) 
into a rectangular matrix with more columns than rows; in other words: 
there will be no square submatrix in A(0) of rank M + 1 if M i N, so 
A(0) shall be rank defficient, and from this follows that A must be singular. 

Emphasis must be placed on the word c omplete when referring to the 
set of "solvable" coefficients: it is possible, by removing some coefficients 
with n < N from the unknowns, to introduce others in their stead with 
n s N, but then the solution will not be a complete set of coefficients. 

(b) 100% aliasing: 

If a term in a Fourier series has a frequency n > N t then it 
will be aliased with lower frequency terms and become impossible to dis¬ 
criminate from them. For most functions of practical interest, the higher 
the frequency, the smaller the term, so the coefficient estimated using the sum in 
(1.21) will be dominated by the lower frequency terms: the estimation error 
is thus likely to exeeed 100%. Estimates above the Nyquist frequency are 
usually regarded as meaningless and the closer a term is to that fre¬ 
quency with increasing n , the less reliance is placed on its estimate. 

In the case of spherical harmonics, expressions (1.16a-b) show clearly 
that the harmonic yf, consists of several Fourier terms of frequencies 
ranging from 0 to n . When n > N , only that part of the Fourier ex¬ 
pansion of above the NyquiSt frequency will become scrambled beyond 
recovery; part of the harmonic is left intact: the low frequency "tail", 
which means that the effect of aliasing on the recovered coefficients does 
not necessarily reach 100% (or even 70%) at the Nyquist frequency, as shown 
in the examples of section 3. 

(c) Orthogonality 

From (1.20a) follows that the matrix of equations (1.19) for the 
Fourier series is orthogonal, so the coefficients estimated according to 
( l 21) are independent from each other. 

Da the case of the Y a ^ orthogonality does not carry over to all the 
sampled harmonics, unless special "quadratures' weights " are introduced 
in (1.5) or (1.6). This lack of orthogonality affects, for instance, the for¬ 
mulas for mean values discussed in section 2. The method of Gaussian 
quadratures is an example of "quadratures with weights" that gives exact 
coefficients when the Nyquist frequency is not exceeded, though it requires 
a special grid where the parallels are situated at the same latitudes as the 
zeroes of P M f 1 (cos9). The use of this method is possible because the product 
[P n ,(cos9) P p « (cos 9)] is a polynomial in cos 9 of degree n + p £ 2N t 
the grid, however, is an unusual one. Details of the application of Gaussian 
quadratures to spherical harmonic analysis are given in a report by Payne 
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(1971). Other examples of methods that recover the coefficients exactly 
when the data is noisless and N, ax < N arc: least squares collocation, 
least squares adjustment, and the algorithm developed by Rice and Burrows 
from expressions (1.17) and (1.18 a-b). 

In general, not all the discretised harmonics retain the orthogonality 
property, and the estimated coefficients are affected by the values of many 
of the other C D ,, in addition to those whose degrees exceed the Nyquist 
frequency. More about this will be said in section 2, when discussing least 
squares collocation and adjustment. 

Summarising: there are enough differences between the aliasing of 
Fourier series and that of spherical harmonics, in spite of their being so 
closely related, to require a great deal of caution before using the intuition 
gained from one type of analysis when attempting the oiher. For this rea¬ 
son, the expression "aliasing error" will be replaced by the just as appro¬ 
priate "sampling error" , which is perhaps less charged with misleading 
connotations because it has not been applied almost exclusively to Fourier 
series. 

It should be noticed here that Gaposchkin (1980) has published formulas 
for the sampling error on a type of equal area grid (i.e., all blocks have the 
same area as the equatorial ones). His formulas are the equivalent, for such 
a grid, of expression (1.22) for Fourier analysis, but much more complicated; 
they are made tractable numerically by the use of certain recursive expressions 
that he provides, thus showing an interesting new approach to the study of 
the problem. 


1.4 Number of Operations in Analysis and in Synthesis 

A 

Expressions (1.17) and (1.18 a-b) can be used to calculate the C„ B 
once the 2-D Fourier coefficients of f(© lf Xj)> or a^, , have been esti¬ 

mated by means of the 2-D Discrete Fourier Transform (DFT). Using the 
Fast Fourier Transform (FFT^ algorithm (Cooley and Tuckey, 1965), the 
number of operations required is proportional to the number of data points, 
which is the order of N s , or CKN 3 ) for short. The FFT is discussed further 
in paragraph 1.9. 

Having obtained the al* B , calculation of (1.17) requires N oper¬ 
ations per c“., or NxN 3 = N 3 for all of them. Finding the coefficients 


^Most of these calculations have the form of scalar products 

N-l N _ 

p = £ a k b k t so the basic operation of finding p k = a k b k + p 1 ' -1 (p° -0, p p) 

k =0 

consists of one sum and one product. 
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1*^ p by means of the recurslves (1.18 a-b) adds another 0(N 3 ) operations 
that can be obviated by computing the I* (P once, and then storing them on 
magnetic tape or disk. One way or the other, cHN 2 ) + 0(N 3 ) operations 
are needed altogether, or cKN 3 ) when N is large (say, N a 180). 

As already mentioned, this prodedure was first described by Ricardi 
and Burrows in 1972; more recently (1977) Goldstein developed a very sim¬ 
ilar idea and formulated a similar algorithm for synthesis. Goldstein's 
method uses recursive formulas for the I* tP different from (1.18 a-b). The 
synthesis algorithm also requires 0(N 3 ) operations. 

As explained in paragraphs 1.5 through 1.7, the procedures presented 
there also require C^N 3 ) operations for analysis, and as many for synthesis, 
though they are formally different from Ricardi and Burrows’. In paragraph 
1.8 it will be shown that synthesis requires as many operations as analysis, 
because one is the dual of the other. 

The fact that two rather different approaches (Ricardi and Burrows' 
and the one described in this report) require essentially the same number of 
operations suggests that 't)(N 3 )" might be a property of all analysis and 
synthesis algorithms on regular spherical grids due, somehow, to the nature 
of the sphere itself. This is speculation, of course, but if not, are there other 
ways of partitioning the sphere for which faster methods exist? The author 
has discussed this possibility before (Colombo, 1979, paragraph 4.6). It 
is interesting to notice that in all these procedures the 0(N 3 ) operations are 
those associated with 0 1( or "column operations"; "row operations" are 
only ©(N 2 ). In the case of the Euclidean plane, the 2-D Fourier transform 
requires the same number of operations per row than per column, 0(N), 
thus the total is only qN 3 ), or O(N) times faster than its spherical "counter¬ 
part." 

While not as efficient as the 2-D FFT, the algorithms for the sphere 
considered here can be much faster than the straightforward implementation 
of expressions (1.5), (1.6), (1.2), or (1.8). The latter has been the ap¬ 
proach of many scientists who have developed their own software, but whose 
main interest has generally been far removed from the study of numerical 
techniques. In 1976, while working at the University of New South Wades 
(Australia), the author developed the two algorithms of paragraphs 1.5 and 
1.6, and C. Rizos programmed them. Subsequently they were used at 
Goddard Space Flight Center, in Maryland. To everybody's surprise, Rizo's 
programs turned out to be more than 100 times faster than those in use at the 
time, when run under the same conditions. More recently, this author has 
written the subroutines HARM IN and SSYNTH described in appendix B. 

SSYNTH has been used, after the fashion of the numerical experiments des¬ 
cribed in Section 2, to generate 64000 1° x 1° mean values (simulated aver¬ 
aged gravity anomalies), each the sum of the 90000 terms of an expansion 


complete to degree and order 300. This took less than 50 central processor 
unit seconds in the AMDHAL 470 V/6-n owned by The Ohio State University 
(OSU). All calculations were in double precision (32 bits words), using 
FORTRAN H EXTENDED. Precomputed values of the Legendre functions 
(or their integrals) were read from tape, and all operations involving trigo¬ 
nometric functions were carried out by a Fast Fourier Transform subroutine. 
These two characteristics, plus a generally tighter coding, are the reasons 
for the greater speed of this program, compared to the older versions men¬ 
tioned before. 

IP all these methods all operations along a given row or parallel (con¬ 
stant 9 t ) are independent from those for any other row, so a parallel processing 
computer with N processors (arithmetic and control units) could analyse or 
synthesize a full grid of N rows as fast as an ordinary computer with one 
central processor can do a single row. This N-fold increase in speed can be 
obtained with the same type of basic hardware (gates, registers) that is 
currently used in conventional "general purpose" main-frame machines. The 
full power of the algorithms presented in this work will be realized when 
computers of parallel structure become more widely available for scientific 
applications than they are today. 


1.5 Algorithm for the Analysis of Point Values 
Expression (1.5) written in full becomes 

A /y N—1 3P^“ 1 < COS t 

C “« [ [ p ». < cos0 i> |sin) mj AX 

1 = 0 J = 0 

which corresponds to the general type 

Att aN “ X aP *-l (cos) 

C a .= L E X a i" E i SinfmjAX fOt.X,) A tJ (1.24) 

1=0 J = 0 i~° 

where x“* could be ^ P BB (cos0t)A tJ as above, or (cos Sj) u>j in 
the case of quadrature w ith we ights u) t , etc. 

To simplify the discussion, the grid is supposed to be equal angular 
and N b »*=N- 1, This and the algorithms that follow can be easily adapted 
for the equal area grids current today. Subroutines HARMIN AND SSYNTH 
(Appendix B) can handle the cases N,» X <N-1 and N,. x >N-1 as well as 
N,„= N-l . 

The equal angular grid is symmetrical with respect to the Equator, 
and assuming that (the same as P B „ (cos 9 1 ) or P tt , (cos ) sin 6 1 ) 

X** = X«-i if a-m is even, and x”' = -X^.i if n-m is odd, one 
can write 


-15- 


c“.= r‘(xv[‘r|sm)^Xf«e 1 .x J )] + 

(■1) xV ^ 21 | sin | m j^^-E(0N-i-t» X a ) \ 


(1.25) 


This formula suggests the following procedure in two nested loops: 
START : set i = -1, ~ 0 for all 0 z n, m £ N ; 

Outer loop : 

(a) increment i by 1 unless i = - 1, in which case STOP 


Inner loop : 

(b) compute all 



(c) find all 


iOt(i) _ 


'a« 


»■ 


+ K 


b 1 . 


+ (-D 


C 1 " 1 1 


v b * 

Xi 


(1.26) 

v 

for 0 s n, m « N (where "(-l) 8-1 | J " merely indicates that is to 

be added or substracted according to the parity of (n-m) ); GO BACK TO (a). 

At the end of the outer loop, = C^, . 


The a 1 , and b 1 , in (1.26) can be computed by taking the Fourier 
transform along row i and row N-l-i of the values of f( 6, X). This in¬ 
volves O(N) operations. There art half as many x“ * as there¬ 

fore it takes 2(N + l) a products, and just as many sums, to form the 

for the pair of rows 1 and N-l-i. Consequently, there are 0(N) + 
0(N*) operations per pair of row^ or 0(N®) + OfN 3 ) for the grid as a whole. 
This is the same as with the Rlcardl and Burrows' algorithm, quickly ap¬ 
proaching Q(N 3 ) as N increases. 


Subroutine HARMEN (Appendix B) Implements this technique. 
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1.6 Algorithms for the Synthesis of Point Values 
Expression (1.8) is of the general form 

N-l a 

f(9i.X 3 ) = [ [ x" [C u cosmXj + S B , sinmXj] (1-27) 

a = 0 ®=0 

where Xi can be, for instance, Kp Bt (cos0 t ), etc., with K being a proportionality 
constant. Rearranging terms and considering the parity of x?" leads to 


\UB it \ s ) _ y /r f 1 ixv 
[f (Wi. X,) “ ,to\Lk |(-D tt -xV 


C M cos m j AX + 


N-l 


[£ l(-s“xv!M 8i ° mi4x ) 

which suggests a procedure in two nested loops: 

START : set i = -1 
Outer loop: 

(a) increment l by 1, unless l = £n - 1, in which case STOP 


(1.28) 


Inner loop : 

(b) compute all 


for 0 s m £ N; 
(c) find all 


1 i N-l I — 

a. = r 1 c n , 

si L ' |s„ 


N-V-l 




N-l 


= Z *"* (-1) 




C.. 


f(Si > Xj) 

» Xj 


N-l 

= Kl| 

* =0 


O. 1 


cos mjA X + 


I2L. 

Pa 


slum jAX 


(1.29) 


for 0 s ] s 2N-1 (where (-1)“"* {—} means the same as in the previous 
paragraph); GO BACK TO (a) . 


At the end of the outer loop all f ( 0,, X 4 ) in the grid are known. 
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Expression (1. 29) is computed by applying the EFT to the 2N ct \, (per row) 
and the 2N 0„~~, taking O(N) operations. The first part of the inner 

loop (forming the , jS, 1 and a",' 1 ' 1 , P!T^) involves C^N 3 ) operations; 
for all £n pairs of rows the total is 0(N ) + 0(N 3 ), and this tends to 0(N‘) 
as N increases. 

Subroutine SSYNTH (Appendix B) implements this technique. 


1.7 Algorithms for the Analysis and Synthesis of Area Means 


(I) Analysis 


Rearranging (1.7) 


% *>■ [&>•.♦ <-« 


N-Ul . 

) 


B(m) 

A(m)] 


b» + ('I)"" 


(1. 30) 


were a., b, are the same as in paragraph (1.5). A procedure similar to 
that for the analysis of point values can be obtained directly by replacing the 
bracket in (1.26) with that in (1.30), and then proceeding as in the algorithm 
for point values. The total number of operations is, once more, OfN 3 ) + 0(N 3 ), 
or Q(N 3 ) for large N. Subroutine HARMIN also implements this algorithm. 


(IT) Synthesis 


6 i+A6 


Truncating the series in equation (1.2), replacing J & P n! > (cos 9) sin0d9 


with Xi’, rearranging terms,considering the parity of x”** and using a, , )9 
as defined in paragraph (1.6), leads to the expression 


Iu 

[i“i- 


- t ! 

- 1 

i 

a a 

! 

I s , 

1 1 

1 

N1-1 
Oi B 

1 A(m) - s ~ 

1 1 

t B(m) 

J 

■ = 0 

L ' 

1 




B(m) 


+ 

$ • 


(1.31) 


A(m) sinmjAX 


Thealgorithm for the synthesis of the F tJ is a direct extension of that for 
point values. The number of operations, once more, is 0(N 2 ) + 0(K 3 ) for 
large N . Subroutine SSYNTH implements this algorithm as well. 


1. 8 Duality between Analysis and Synthesis 

Pairs of direct and inverse linear transforms, such as Fourier trans¬ 
forms, possess dual characteristics: certain words and mathematical 
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expressions can be arranged in pairs (a,b) such that, if every "a" is replaced 
by its "b" in any statement or equation valid for a function f , then the mod¬ 
ified statement is valid for the transform F c£ f, and viceversa. 


Analysis and synthesis of spherical harmonics are reciprocal linear 
transformations of data into coefficients and of coefficients into "data" closely 
akin to Fourier transforms, so they can be expected to exhibit some dual prop¬ 
erties. Comparing the formulas and algorithms in paragraphs (1.5), (1.6) 
and (1.7) shows many similarities, among them the number of operations. 

This can be understood as being a consequency of duality. To make this point 
clear, consider the following pairs: 


( 

.Xj) 


( 

A 


f( ®i» * j) 1 )(and similarly with 

ft: 

V 

f(®v_i_t X j) 

■ 



’ 

f(0N-l-l» ^j))/ 

Th-i-i 


(v a * /. 

\X i r 


l) n_ * , sinmjAXi ; (i, m) 


(3.a); (X?. 


cos mjAX) 



From these we can derive the following pairs: 


a* 


a' 




N-l-i 


fe-! • 


Pm 


A(m) 

ffi(m) 


a 1 . 


B(m) 

A(m) 


A 
■ pT 


(a 1 , , A(m)); (b 1 , , B(m)) 5 (aT 1_t ,-B(m)) i (b* 1 1 , A(m)); etc., 

and, in conclusion: 

("ANALYSIS", "SYNTHESIS") 


Each one of the analysis algorithms becomes its synthesis counter¬ 
part by a simple replacement of terms. Once an algorithm for analysis 
( synthesis) is defined, the corresponding algorithm for synthesis (analysis) 
follows. For instance, one can easily apply the principle of duality to the 
Rlcardi and Burrows’ method of paragraph (1. 4) to obtain a synthesis technique. 


1.9 Usefulness of the Fast Fourier Transform Method 


The excellent book by Brigham (1974) gives a thorough presentation 
of the l-D discrete Fourier transform and its applications, and explains 
in detail the method known as FFT for computing such transform. The 
Fourier transform in 2 and higher dimensions chn be found simply as follows: 




first, get the 1-D transform of each row, then that of each column of the mod¬ 
ified array .... etc., until all dimensions have been exhausted in this way. 
Understanding the workings of the l-D transform is enough to understand 
those of the N-dimensional transform as well. The FFT requires 0(number of 
points) operations for each row along the nth dimension, so the number for 
all points in a regular, euclidean array is always of the order of the total 
number of points in that array. 

Before the mid-sixties _when the FFT came along _ the best techniques 
available for the analysis of data on regular arrays required 0(number of 
points) 2 operations. The increase in speed of o(aumber of points) brought about 
a true revolution in data processing: work that had been long regarded as 
impossible became feasible overnight, the field of industrial and scientific 
applications for numerical Fourier transforms expanded tremendously; 
the impact in areas as diverse as cristalography and communications engi¬ 
neering was remarkable. 

Having mentioned the positive side of the FFT, which is used in the al¬ 
gorithms described so far (at least in principle) and in the programs HARMIN 
and SSYNTH, it is only proper to say something about its alternatives. 

The FFT calculates a!12N Fourier coefficients a,, b, very efficiently, 
but takes just about as many operations to get only a few coefficients as it 
takes to get all: forN,,* small compared to N there may be a real dis¬ 
advantage in using the FFT. The FFT is most efficient when the grid is 
such that N is an integer power of 2. The grids used in geodesy are usually based 
on the division of the circle In 360°, and many on the sexagesimal division 
of 1° as well, hi all of these N contains factors other than 2, so a less ef¬ 
ficient version of the FFT, known as the mix-radix FFT (Singleton's algorithm) 
must be used. 

Finally, the mix-radix algorithm is rather convoluted, so it is best 
to take ready available subroutines from software libraries (as it is done in 
HARMIN and SSYNTH) rather than to incorporate the FFT "on line" in the 
program one is writing. This means that the program is going to be less 
self-contained. 

The "pre-FFT" methods can be more efficient than the FFT when 
N,,«« N; they are also very easy to program. For the sake of completeness, 
the outline of a method this author has used quite often will be given here. 

Consider the trigonometric relationships 
cos(aj9) f 2cos,8 cos(a-1) 0-cos(ot-2)0 (1.32) 

sin(a0) = 2cos (3 sln(a- l) 0- sin(a - 2) 0 (1.33) 
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with a, j8 real. If all the values of the trigonometric functions could be 
obtained with one operation oaelu the number of operations involved in 


find ing 


a N-l 

'So 


cos 

sin 


mj AX (X s ) f( 0j, Xj) woulc3 be 2 ( 2N )» or 4 ^' 2 for 


all 0 < m ^ N. In other words: 0(N 2 ), as expected. This is precisely what 
can be done using (1.32),(1.33) as recursive expressions for cosm(jAX) and 
sinm(jA X) with m integer and 0 £ j s 2N - 1. The values of 
cosm(-AX) = cos mi X and sinm(-AX) = -sin mA X , are needed to start the 
recursion; they can be calculated with standard trigonometric subroutines. 

The use of such subroutines increases the number of operations slightly over 
4N 2 , but if N is large enough, this is negligible. 


With all calculations carried out in double precision, this method 
gives values of cosine and sine that coincide to better than 5 significant 
figures with those provided by the standard FORTRAN functions, when 
N is as large as 1800 (0.1° x 0.1° grid). By taking advantage of half-wave 
symmetries in the sine and cosine, and by ingenuous programming, the num¬ 
ber of operations can be reduced by a factor of 4 or more rather easily. 


1.10 Functions Harmonic in Space and their Gradients 


If f(0, X.r) satisfies Laplace's equation V 2 f = 0 in the space outside 
a sphere of radius a , then it can be represented, in that space, by the 
solid spherical harmonic expansion 


f(0,X,r)« tit Si C »■ Y “ 

n ~ O a — o uT= 0 1 


(1.34) 


If we consider, at a point P ^(0,X,r) in space, the local triad r, h, t 
oriented downwards to the origin, West to East along the tangent to the local 
parallel, and North-South along the local meridian, the components of the 
gradient of f(0, X.r) along this three axes are 


M (9 

3r 


3h 


n = o B = 0 

co n 


n = o » = 0 

m C na sin m X] 

oo n 

|f(0,x,r)= £ £ 

at n=0 "=0 


ll 

o 

a n , , —a 

pr+s (n +1) C n9 Y n „ 

(0, A) 

(1.35) 

a n 

yTC+a 

P nB (cos 9) cosec (8) 

lmS BB cosmX - 




(1.36) 

a n 

]T» +S 

dp n» (cos0) IC„, cosmX+5' 00 
d 0 

sinmX] 

(1.37) 


1 


^ereonCoperatiorf,' as mentioned in paragraph (1.4), consists of one 
sum and one product. 
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Expressions (1.34), (1.35), (1.36), and (1.37) appear often in the discussion 
of geodetic problems, and their calculation from a set of coefficients is an 
important problem. If all the values of any of therewith r constant^nd 
on a regular grid) arc required, the methods discussed so far can be used 
after a few minor additions. 


The synthesis algorithms can be thought of, for the purpose of this 
discussion, as "black boxes" with the coefficients C bB > the Legendre func¬ 
tions, N„ ax , and N as inputs, and all the 2N 2 values of f(8, X) on the 
corresponding regular grid as outputs. To compute the expressions given 
above, only the part of the input consisting in the coefficients and/or the 
Legendre functions has to be modified before they enter the "box", which 
r^pains untouched, ^or instance: to compute (1. 34) one should replace 
•C,with " pj +2 C 0 „ " in the "input"; the others are equally obvious and 
will not be explained further. The following recursive formulas can be used 
to obtain the derivatives of the Legendre functions: 

dP„ = (sin0) -1 { n P nB (cos0) cos0 - 1 ^ P n _ x m(cos 0)f 

d“0 L < 2n - V J ( i.38a) 

f-f nn = l sin0 -- 1 + cos0 P *- J - 1 } (1 * 38b) 

with the starting value 

%°° -° 

These recursives follow from the unnormalized formula 


(cos 2 0- 1) d_P^ B (gos 0) = n cos0 P BB (cos0) - (n + m) P^^cos©) 

(N.N. Lebedev, "Special Functions", Dover, 1972, Ch. 7, equation 7.12.16), and from 

i 

P = [ il 1 ) 1 2 sin0 Pn-i n-1 (see paragraph (4.4)) 

L (2n - 1) J 


and 


- [ - <n /; !! - , P P.. ; - /STTi p. 

L 2{2n (- 1) (n - m) ! J 


2(2n (- 1) (n - m) 

The complete recui'sive expressions for the P B a are given in paragraph (4.4). 


The expansion for the area means defined by (1.2) can be differentiated 
terra by term because it converges uniformly. The expressions for area 
means gradients, equivalent to those given here for point values^ arc immed¬ 
iate. They can be computed after simple modifications to the C B0 and/or 
the xV* :md usin S the same programs for computing the area means. 

Subroutine "LEGFDN", listed on Appendix B, can compute both the nor¬ 
malized P nB (cos 0) and their derivatives dP B „(cos0). 
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2. Error Measure and Optimal Quadrature Formulas 

This section introduces a criterion for quantifying the errors of num¬ 
erical quadratures formulas that is based on the statistical properties 
of the data. Three qualities are highly desirable in an error measure: (a) 
it should be easy to determine; (b) it should be mathematically tractable; 

(c) it should provide a good idea of the likely size of the actual errors. Point 
(a) is taken into account by choosing a quadratic measure, because the numerical 
formulas are linear estimators of the C®,, and the linear, quadratic esti¬ 
mation problem is fairly simple, with its mathematical side very well under¬ 
stood and developed today, which takes care also of (b). Regarding (c), the 
reader will have to wait till section three, where certain evidence, obtained 
from numerical experiments, supports the assertion that, though statistical 
in nature, the error measure adopted represents the actual errors very 
closely. 

Having defined the error measure, the notion of optimal or best formula 
according to such measure is investigated, leading to the application of least 
squares collocation and least squares adjustment to spherical harmonic analysis. 


2.1 The Isotropic Covariance 

The Isotropic covariance (expression (1.10)) between two functions u(6,Ai, 
v(9 ,A) on the unit sphere, both expandlble in spherical harmonic series 


00 n 1 

u(0. A) = K u £ Z Z 

aft 

^(d,X) ; 

n = 0 »=o 0 P= 0 



v(6. A) = K «Z Z L 

^ n* 

X) 

a = 0 ,= 0 0=0 




can be formally defined as follows 

cov(u(P), v(Q)) = Mju(P) v(Q)j (2.1) 

where M { } Is the isotropic averaging operator and P and Q are two points 
on the sphere separated by the spherical distance 0 t9 . The operator m{ } 
symbolizes the average of its argument (in the present case the product 
u(P) v(Q)) over all rotations of the sphere. This can be visualized if one 
thinks of the points P and Q as given in a fixed system of coordin.tes, while 
the sphere, on which u and v are defined,rotates In all possible ways. After all the 
(infinitely many) possible rotations, the average product u(p) v(Q) will be iden¬ 
tical to cov(u(P), v(Q)). This kind of covariance, though purely geometrical, 
resembles closely that of stochastic processes such as time series. 





The importance of the Isotropic operator and the isotropic covariance 
function in spherical harmonic analysis stems from the fact that the latter can 
be described as the estimation of certain parameters of a function f (S, X), 
the Cj, , from data sampled on a sphere. From paragraph (2.7) on, this 
report deals with optimal estimators for the C *, based on the theory of least 
squares collocation. Such optimal estimators minimize a quadratic measure of 
the error that is defined in terms of the operator M { [ , this measure 
being introduced in paragraph (2.4). 

The idea of least squares collocation is related to the basic principles 
of such linear, minimum variance estimators for time series as the Wiener 
and Kalman filters, which have found wide application in the physical sciences 
and in engineering over the last thirty years, and have been generalized 
to deal with both continuous and discrete time processes, and also "processes" 
in more than one dimension, such as are found in pattern recognition and in 
digital image enhancement. Two-dimensional Wiener filtering, of which the 
reader can find several fine descriptions in the special issue of the Pro¬ 
ceedings of the IEEE, Vol 65, No. 6, 1977, is also applicable to "flat-Earth" 
geodetic calculations; least squares collocation can be regarded as the ex¬ 
tension of this type of filtering to calculations on the sphere. Isotropic av¬ 
erage operators are not the only ones that could be used in the "statistical" 
approach, though they are probably the easiest to work with and, perhaps, the 
best for the sort of application considered here. For a description of other 
likely operators, the reader is referred to the paper by Rummel and Schwarz 
(1978). Probably the most didatic introduction to the method of collocation 
remains Heiskanen and Moritz, (Ch. 7, 1967). 

Reasoning as in Heiskanen and Moritz (ibid), one can show that 

CO 

CQV(U(P), v(Q) ) » Z C “n ,V P»(COS^pq) 
assO 

which is, in fact, expression (1.10*) _ the definition of the isotropic covariance 
given in section 1 without any reference to M { } . Similarly, 

00 

COV(U(P), U(Q))= Z <*n Pa (COS 0 pq ) 
a= o 

usually known as "the covariance of u " (expression 1.10), while (1.10*) 
represents the "covariance between u and v ", or "the crosscovariance 
of u and v ". The one-to-one relationship between covariance and power 
spectrum (or crosscovariance and crosspectrum) should be clear 
from these expressions. 

—CK 

To apply the notions introduced above to the C a ,, it is necessary to 
think of them as functions rather than fixed values. This is possible if one 
considers changes in the coordinates 8, X brought about by rotation. Each 
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such change results In different coefficients, though the function they describe 
is the same, only rotated with respect to the old system. The new system 
can be related to the first by three angles: the coordinates 9, X of the shifted 
north pole, and the azimuth A of the zero meridian. Therefore, the C® t 
are functions of 9, X, A and this is enough to define the average over all 
rotations of the product of each C®, by itself or by another function, in a 
meaningful way. Two important properties of spherical harmonics are: 


(A) 


I A — a — g 

M £ c„ + s„ 

( * =0 


M 



a a 2 


( 2 . 2 ) 


i.e., the power spectrum is invariant with respect to rotations. This 
follows from the plain fact that the integral f f a (0, X) da is invariant 

over rotations, and from Parseval's identity (1.12); (2.2) Implies that the 
isotropic covariance function (1.10) is likewise invariant. 


(B) 


M 


_a _ p 

Yb.(p) Y,,(q> 


= 0 if 


a t 9 
n t p 
m 4 q 


for all P and Q 


(2.3) 


i.e., the orthogonality properties of spherical harmonics with respect to 
integration on the sphere are also true with respect to averaging over ro¬ 
tations. 

As a consequence of (2.2) and (2.3) above, the following relationships 


are also true: 



M|c„. 

(2.4-a) 

sy a 



M c!,| - o if j 

n 4 p 

(2.4-b) 

1 

m i q 


Mjcf. ( * 0 if n 4 0 

(2.5) 

M|c“. f(0,X)} * ** (0 ’ X) 

(2.6) 

M I C a ■ 1 u| (in + l)h 

ITAi Y®(9,X)da 

(2.7) 
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For the derivation of (2. G) the reader can see (Ruramel, 197G), and (Sjoberg, 
1978) for (2.7). As for (2.4a-b) and (2.5), the proof is given now. 


According to (1.4): 

16Jr 8 m{(c“,) 2 }= m {/^(0^) f(e;x<) dajf (e;x?) f<e;V) da*} = 

* M {^ Y?. (0,X)f(0, X) . f ( 6', X') Y^j, (0 \ X') dff do’ = 

f a V n a B (0,X) y“ (0\X«) M{f(U)f(S',X')} dada’ = 

fa fa' Y® (6. X) y“ ( 0',X') cov(f(P), f(Q))dada* 
where P=(0, X) and Q = (0',X'). According to (1.10): 

GO 

167T 2 m{(C^ b ) 2 }= fa X“<®.X)da f a , y“ (0',X') I o a 2 P n (cos,k,)do’ = 

a ' 2 ^2 

=4TT Jo Yn* ( 0 > X) Y® (0 , x ) da = 2 ^— 16U 2 because of ( L. 14) 

Similarly, 

r _tt —9 i a,? i r — a —b ® ^ ^ 

M[c n .C P ,}= gn+T * "^7 J<J Y an (0,X) Y p p q (0,X)da =0 if n/p 

m^q 

Finally, recalling that Y^ o (0, X) = P 00 (cos 0) = 1 for all -tt £ 0 £ n and all X, 

m{c*}= pi / a Yf,(0,X) M{f ( 0,X)}da=M {f(0,X)} Y® (0,X) da 

= M.{ f < 6 ’X)}^ y“( 0,X) Yg o (0,X)da =0 IfnM, 

which completes the proof. 

2. 2 Some Additional Notation 

So far, data points on the sphere have been identified by the subscripts 
i and j . Alternatively, they could bo arranged according to a siigle subscript 
k = 2Ni + j (where N is the number of parallels in a grid and 2N the num¬ 
ber of meridians), so the points in the "O' 1 row, ordered by increasing j , are 
followed by those in the "l" row, in the same order, etc., the last element in 
the "N - 1" row closing the sequence. Based on this convention, the set of all 
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values of f (0 1 ,Xj) (or f i} ) can be arranged in N p - vector form according to k 


I z o z r ■ • z i 


z Nrll 


where 


z k ~ or z k - T t j 


with k = 2Ni + j 


and where N p is the number of data points, or 2N 3 for equal angular grids. 
In a similar way, the coefficients C „ B can be ordered according to a single 
subscript p = n 2 + an + m + 1 (with the understanding that the meaningless 
5 n0 are not included ) defining the following N a - vector: 


[ C 0 c x . . . Cp . . . C N(1 ] T 


Co ~ C. 


p = n 2 + an + m + l 


where N 0 = (N b0X + l) 2 . Using this notation, expression (1. 27) for point 
data quadratures can be written 


= Inn £ 


( 2 . 10 ) 


where the estimator vector f aB , of dimension N p , has elements of the form 


v n« {cos'! 
l sin/ 


mj AX 


under the convention given above relating k > I and j.. 

a a a 

Grouping all the estimates C BB in a vector c ordered in the same way 
as c , the relationship between the and z can be written, in matrix 
form, 


c = F z 


( 2 . 11 ) 


where F is the estimator matrix implied by(2.10). It is a N p xN 0 matrix 
(where N p is the number of data points in the grid), each row bein^formed by 
the coefficients of the quadrature formula for the corresponding C nn . Such 
row is also the transpose of the estimator vector of this , designated _f^ 0 

in (2.10). 


In the same way as the covariance function between scalars, the cov- 
iancc between vector functions can be defined in terms of m| j-: 


r n 

M ' L Z £ j = C lt 


( 2 . 12 ) 


where C tz is the covariance matrix of of dimension N p xN, , This 
matrix if a function of the relative positions of the points in the grid on which 
?_ has been determined, in the same way as the scalar covariance depends only 
on the distance between two points. The elements of C rI arc 







C «\ = M ;z f z.} = M;f(P r ) f(P.)}, 

i.e., the values of the scalar covariance corresponding to pairs ef points 
in the grid. 


In the same way 

f 

M\cc J- = C (2.13) 


is a N a x N a diagonal matrix according to (2.4a-b). C is the covariance 
matrix of the coeffic ients. Similarly, the covariance between c and z is 


M 



» C 


0 2 



(2.14) 


where C tt is a N a x N p matrix, the elements of which are 

c «| = M i_c p z,} = M^C tt> f(Si,Xj)j' 

where the right hand side is given by (2.6). 

Finally, when estimating the cf B , not from samples of f (9,X), but 
from measurements corrupted by noise 


mtft.Xj) = f(«t,Xj) + n n (2.15) 

the measurement errors can be grouped in a N p - vector n with the same 
ordering as z , and the sum of both will be, then, the N p - vector of ob- 
served values 


m = z + n 


(2.16) 


The measurement errors are values that occur in time, as successive 
observations are carried out: they constitute a time series. The average 
operator appropriate to them is the usual statistical expectation operator 
E ^ j . The measurements are suppossed to be unbiased , so E ■ n k j- = 0 
for all k . The covariance implied by this operator, is theusual statistical 
covariance: E|n k j « cfc 3 * , and E |n k n r | = . This can be 

generalized for the noise vector n : 


E{nn T } “ D 

where D is a N, xN, matrix of elements 
d *r - E{n k n,} = <r* r 


(2.17) 


(2.18) 


Both C , , and D have in common a very important properly: 
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x T C,,x ^ 01 If ?c T ^ 0 , x any N p vector, 
x T D x 2 0 J 

i.e., they are always positive matrices, moreover. In all the cases considered 
here, at least D is positive definite ? 

x T D x > 0 for all x 

Positiveness can be Inferred readily from the definitions of C tl and D : 
i.e., x T D x = x T E J L nn T }x = e{x T o t x} = E {h 2 } *0 

(where h = x T n), and similarly for x T C tt x(with j) . 

2.3 Estimation Errors, Sampling Errors, and Propagated Noise 
A linear estimator is of the general form 
£ = F m 

where m is the vector of measurements defined by (2.16), and £ is the vector 
of estimates, made up in our case of the (5 a °^. According to (2.11) and (2.16) 

£ s c = F (Js + n) 

hi general, the estimates will not be exactly equal to that which is estimated, 
the difference being the estimation error. In matrix notation 

e = c - c = (c - F z) - (F n) (2.19) 

e being the estimation error vector. The two terms in the expression above 
can be defined as the components of this error: 

e, = £ - F z 

which is the estimation error in the case of noisless (perfect) data; and 

2.V = En 

which is the error due to the noise, or propagated noise . 

The error e Ip may be due to a number of reasons. 1 If it is zero for 


1 Using the relationship p = n a + an + m + l of paragraph (2.2), e, p stands 
for the sampling error in 6* m , and e^p for the propagated noise. 


some estimator, then its presence in other estimators could be blamed on 
them being somewhat inadequate. For Instance, if the estimator was chosen 
by taking the elements of F from a set of random numbers, then the estima¬ 
tion error is likely to be always high, as the estimator has nothing to do with 
the actual problem. In particular, the addition of extra measurements to 
the vector m is not going to bring any general improvement on the estimates. 

On the other hand, if attention is paid to the nature of the problem when selecting 
F, one would expect the error to decrease as more data is introduced. If, 
as the number N, of samples in m tends to infinity, e*j, tends to zero, 
one could say that the error-is due to the incomplete sampling of the signal 
f ( 9, X), and call it the sampling error. This is precisely the case with 
any of the quadrature formulas to be studied here, all of which can be written 
formally as linear estimators according to (2.11), and for all of which the 
error e, 9 vanishes as the number of samples tends to infinity, because the 
sums become identical with the Integrals defined by (1.4). In this sense it is 
quite suitable to call e, p the sampling error, as in paragraph (1.3). 


2.4 The Quadratic Error Measure 

The overall error measure will be defined here as the sum of two quad¬ 
ratic terms: onefor the propagated noise, the other for the sampling error. 


(a) Propagated Noise Measure 


This measure is the same as in least squares adjustment, i.e., the 
variance of the error defined in terms of the usual statistical expectation 
operator 

= E{(f“ T n) 8 } (2.20) 

according to (2.10). This variance represents the scatter in the value of 
<5® due to the uncertainty in the values of the data. In matrix form 

Etj = E {e^e^} = E {Fnn 1 F T } = F e{u T } F T = FD F T (2.21) 

where Efj is a N„ x N, matrix, while D was already presented in paragraph 

( 2 . 2 ). 


In the special case where the measurement errors are uncorrelated, 
D is diagonal, and (2.20) becomes 


a T , f.OlT t _Q£ "* ,ar a 

° rf ~ E il“« ° 2. in« J' ~ lr.a D_f n , = 


V C COS a \ , 3 ^ 

2. (x“") l Uin a i m ^ XE i“w/ 

1 =0 3 = 0 

which is the usual formula for propagating the covariance of the noise. 


( 2 . 22 ) 


(b) Sampling Error Measure 


This measure is defined in terms of the isotropic averaging operator of 
paragraph (2.1) 

a. 2 „f= m{(c p -f“ Z ) 2 } s M {(c“ ) a } (2.23) 

or, in matrix form 

E s = M-fe.elj- = M {(c - Fz)(c - Fz) T } (2> 24) 

= C - 2 C e , F t + F C zz F t 

where E, is a N a x N e matrix, and C , C 9z , C zz were introduced in 
paragraph (2.2). 

(e) Total Error Measure 


The total measure is the sum of (a) and (b) 

_2(V _ J2 rj - . „2Q; 

or,« cr__ + cr, .. 

r. » f] 8 r - s 

or, in matrix form, 


E r = E s + E^ = C -2C e2 F T + FC n F T + FD F T = 
C - 2C az F T + F(C„ + D) F t 


(2. 25) 


( 2. 2G) 


where E T is the N p x N p error matrix associated with F and with the 
covariances that define C , C cz and (C„ + D). Expression (2.2G) is a 
special case of the formula for "E, 6 " in least squares collocation (for instance, 
Moritz (1978), Ch. 3, eqn, (3.20)); moreover, it belongs to a family of 
formulas also found in the minimum variance estimation and filtering of time 
series and of processes sampled on the euclidean plane. 

The total measure has been chosen simply as the sum of 
by making the basic assumption that the sampling error and the propagated 
noise are due to completely independent causes. The first depends on the values 
f(0i»X j )> while the second depends on the measurement errors of instruments 
that, at least ideally, operate with accuracies unaffected by the quantities 
measured, or in such way that any interactions can be eliminated by simple 
corrections. 


The columns of F are defined by the quadrature formula used, and 
such formulas either satisfy, or tend to satisfy, orthogonality conditions (para¬ 
graph (1. 3 (c)). For this reason, provided that C IZ and D belong to the type 
to be described in paragraph (2.9), matrices E^, E r . , and, thus, E T , arc cither 
diagonal or diagonal dominant, and in the latter case tend to become diagonal as 
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the sampling intervals decrease, or N p -* “. For this reason the correlations 
among the eri’ors for individual coefficients are, or "tend to be", very small. 

The diagonal elements of the error matrices E s , E„ and Er are the 
variances of the errors in the respective coefficients, as defined by (2.20), 
(2.23), and (2.25), respectively. 


2.5 The Meaning of the Error Measure 

The treatment of the propagated noise is the same as in least squares 

adjustment, so this part of the error measure should be easily under¬ 
stood. The sampling error measure, on the other hand, is a geomet¬ 
rical measure; M { } belongs, as a concept, in the field of integral geometry, 
or the study of "geometric probabilities". This is a branch of mathematics 
closely related to integration and to measure theory, and also to statistical 
mechanics. In geodesy, this type of idea is relatively new (Kaula, 1959), 

Moritz (1965), but it has been used already extensively enough to show its 
considerable worth. 

From expressions (1.10) and (1.11), the covariance and the power spec¬ 
trum are functions of each other. Since either of them, and the sampling 
grid, define matrices C , C 8Z and C ZI in expression (2. 24), it follows that 
a statement on af aa is, somehow, also a statement on the performance ;>f F 
for all the functions that have the same power spectrum that determines the 
diagonal elements of C . To put this more precisely, consider a function f x (0,>.) 
having the given power spectrum. If C n °l were estimated for f x and 
also, at least ideally, for all its rotations, then the mean square of the sampl¬ 
ing error e Jp in C a ® for all this functions would be, by definition of M J , 
the measure a aa ^ . If a second function f a (perhaps not a rotation of f t ) and 
all its rotations were then analysed In the same way, the average of ef D for 

? /y ' 

all these functions would be, once more, a, nB , as long as f 3 has the same 
spectrum as f x . Moreover, the average of for fj , f 3 , and their ro¬ 

tations put together , would also be crf aB . In fact, if we had a finite set of 
functions I 1 , f a , . . . f B , with arbitrary n , all with the same power spec¬ 
trum (or covariance), then ef p would average crf a ^ for all the fi and their 
rotations. 

It appears, from the preceeding discussion, that one could take a simple 
step and say is the mean of the sampling error squared of the estima- 

tor <5 b ® = f £ , over all possible functions with the given power spectrum. " 
Unfortunately, as mentioned in the Introduction, the sphere is a rather 
wicked surface. There is a theorem by Lauritzen (1973) that states the im¬ 
possibility of having the same average crf BB for all functions as for evei'y 
function, when the distribution of the ensemble happens to be gaussian. Moritz 
(1978) has endeavoured to show that this is no problem if the ensemble of functions 
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is not gaussian, hut using his conclusions here would force the introduction 
of a rather strange requirement of "non-gaussness" on the ensemble of the 
signals analysed that is best left out, if possible. 

Perhaps there is a way out in going back to the idea of a finite set of func¬ 
tions f t , where the problem does not exist, by saying: 

"the error measure a 3 2 n ^, for a certain estimator and a certain typo of 
signal power spectrum, is the average of the square of the sampling error in 
for all functions with the given power spectrum KVKK TO RJC ANA J/y.SKI) 
with that estimator, and for all their rotations". 

After all, accuracy is what geodesists are always interested in, not 
perfection. 


2.6 Simple Formulas for Area Means 

The numerical studies of section 3 concentrate in area mean type for¬ 
mulas, because area means are preferred for collating information, 
particularly on a global basis, at present. The formulas to be studied here 
and in that section can be divided into "simple" and "optimal". The name 
"simple" is given here to expressions of the type 
Aa n-i am a _ 

C na = |i» [ [ f», / Xi. (^A)dcr ( 2 . 27 ) 

l=0 J= o J 3 


where fi n in a scale factor affecting the nth harmonic as a whole. Ex¬ 
pressions of this type have been developed more or less intutitively, along 
the lines of the following reasoning: 


If the signal were constant on each block, it will equal its mean value 
there, and the coefficients of such a fuention would be precisely 



1 

4 TT 


N-l 

I 

1 = 0 



(M) 


(2. 28) 


according to (1. 4). In general, most signals are not equal to their mean value 
over whole blocks, so the expression would not be exact. In most eases, 
the signal would have fluctalions in each block, and it would be less smooth 
than a function that is constant overreach block, so using the formula above 
with Til as data may result jn,the C® of a smoothed function. As a refinement, 
one could try to do-smooth the C n ^ . If the blocks were circular, the relation¬ 
ship between "true" and "smooth" TT n $ would be 
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(2. 29) 



where /3„ is known as the Pellinen smoothing factor of degree n 
The relationship between /3 n and the radius of the circular blocks is given in 
paragraph (4.3). For small blocks, experience shows that there is little 
difference between the area means of geodetic data on circular or on square 
blocks, so the error is small if one assumes that they are the same; in 
such case the modified expression 

N-l 3N-1 




4nj8 n 


i-o 


I fn 

i~ o 


/ V-“ 
•>Oi j 


Y na <e,X)dcr = e 


o! 


(2.30) 


could be used; in practice, this is only an approximation,though a good one, 
as showed by Katsambalos (1979), who tested this expression extensively. 


In addition to (2.28) and (2.30), Lowes (1978) has proposed using 


A 

C 


a 

na 


4ttV 


N-l 8N-1_ » n 

l l U, Y nB (0,X)da 

1 = 0 3 = 0 


(2.31) 


to estimate the harmonic coefficients. All these expressions have the prop¬ 
erty that, because - 1 , and } }^nfa iS ( e »^) da f(0,X) Y a (0,X)da 

as - 0 (or N p - 00 ), it is true that the error e tP = c® - - 0 with 

N P -* <» ; in other words: e ap is properly called a sampling error in the sense 
given to this term in paragraph (2.3). 


Comparing (2.28), (2.30), and (2.31) it is easy to see that they all belong 
to a class of expressionsof the form (2.27), with /j n = , /i n = , and 

Md = JTFg'a , respectively. 


The scaling factor /j. u can also be regarded as a de-smoothing factor, if one 
wishes to retain the intuitive meaning of these formulas. In the notation of para¬ 
graph (2.2), these expressions can be written, according to (2.27), as 




nm 


a T 

Mnd2.no) {l 


with 




a 

no 



(2. 32) 


Replacing (2.32) in the definition of the sampling error measure, (2.23), 
and adding with respect to m and a to obtain the total error in the nth harmonic; 

£ Y_ ft» = °’ 2_ [ 2 Y. — »"0f * ilo*] Mn + 
n=o,=o La=o , = o J 

r £ [ <h%J J (C tt ) h?l /4 2 (2.33) 

jp£— 0 B -0 J 
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r 


(where £},*&, i is a row of C ci ). This is the sum of certain diagonal ele¬ 
ments of E s , according to (2.24), when the estimator has the form (2.27). 
Clearly, (2. 33) is a quadratic function of the scalar , and as such it can 
have either a maximum or a minimum. If C IZ is positive definite, it must be a 
a minimum. Finding the corresponding value of p n is the same as finding 
the formula of type (2.27) that has the smallest sampling error per harmonic 
for signals with the covariance (power spectrum) specified by C zl . In ad¬ 
dition to the sampling error, the measure of the propagated noise can be added 
to obtain 


I I orf- 

Z I 2 

a ^ 0 

?=o, =0 

r ' ” , a t 

. a' 

I Y / (linn) (Czi 

D) h s0 

0 

0 

III 

0 

- 


(2.34) 


This is also quadratic and has a minimum, and finding the optimum p n is 
the subject of the next paragraph. 


2.7 Optimum de-Smoothing Factors 


The coefficients of p n and p n in (2.33; and (2.34) are both real 
scalax’s, and so is the independent term of . The expressions represent 
parabolas, and because both C zt and D are positive, if the further (and 
likely) assumption is made that they are also definite, then b' ! V J > 0, 

Hu? <T « 


and the parabola has a minimum where p n satisfies the condition 

l I otZ “ -I lc_U, h“ * 

« b <x s 


1 £. 

2 Bp, 






a 

Mn 


for the sampling error (2. 33 


i.e., at 


or at 


P a = 


* 

Pn 


^"£113^1 jinn 

l J(h")C„h^ 

\ X £»»«,» ii«a 


(2. 35) 


( 2 . 36 ) 


l U hV(C. + D)h“ 

for the total error (2. 34). 

Expressions (2. 28), (2.30), (2.31), and (2.36) will be studied further, 
by means of computed examples, in section 3. 
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2.8 Least Squares Collocation 

In the notation of Paragraph (2.2), optimizing /i a is the same as ob¬ 
taining the optimal vectors _f of the form 

.a 

for the estimator 

A _ 

£ = Fm 

where m = £ + n , and the (f^) T are the rows of F . If no restriction is 
placed on the form the rows of F can take, then a reasoning similar to 
that in the^areceeding paragraph leads to the best possible linear estimator 
for the C„ . 


Considering the total measure of error for 0 sn sn : 

I N-1 « „ N-l I W-l a n 

»* - E E t < -l<-* HZ 

a=oa=o«=:o a=o a=o ®=o *=o 


+ 


r«. T 


' N-l n 

III (*»>' (C„+ D)f 
a=o u=o«=o 


a 


(2.37) 


it is not difficult to see that, because all of, are non-negative, finding the 
F that minimizes their sum is the same as finding the F that minimizes 
them individually . The sum of the mean squared errors of all coefficients is 
the trace of the error matrix E T of (2.26) : 

- t T t <4® - tr[ E t ] 
ot=o» = ° *=° 

To obtain the condition for a minimum, one must differentiate (2.37) so, 
according to (2.26), 

\ [Et] = -C l t +(C„ + D) F t (2.38) 


as found using well-known matrix analysis formulas. 

From this follows that 

F - C, t <C„ +D)' 1 (2.39) 

is the F that minimizes (2.37), provided that (C„ + D) Is positive definite . 

As already explained, both matrices are always positive and their sum is 
usually definite. The expression for the optimal estimator for is 

= (ff,) T m = (C w + D f 1 m (2.40) 
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where (£“) T 
6 ®. 


is the row of the optimal estimator matrix F corresponding to 


The use of expression (2.40) is, in brief, least squares collocation 
applied to spherical harmonic analysis. 

When the optimal F is used, the error matrix becomes, according 
to (2.26) and (2.39), 

Et — C - Co* (Cjj + D) Cox (2.41) 

and the total error measure is the trace of this rratrix: by definition, the 
smallest for all possible F. 

Clearly, whether one is interested in estimating coefficients or in 
determining the likely accuracy of such estimates,, using expressions (2.40) 
or (2.41) require a knowledge of either (C„ + Df 1 (inversion) or, at least, 
of C 3 j(C„ + Df 1 (solution). Because the matrix (C* r + D) has dimension 
xN s , obtaining either requite s, by usual linear algebra methods, O(Np) 

(or 0 (N 6 )), operations. In the case of a 1° x l°grid, N 5 = 64800, so, at some 
200000 products and sums (double precision) per second, a modem computer 
like the one at OSU would need about one century to obtain all 6?, to degree 
and order 180 from data on such a grid. Fortunately, as explained in para¬ 
graph (2.9), if the covariance functions of signal and noise both satisfy certain 
conditions, and if AX is constant for the whole grid, then both Cxx and D 
(consequently their sum) can be inverted in much fewer operations than by con¬ 
ventional methods, because they possess a particularly strong structure. More¬ 
over, the optimal estimator C ® m turns out to be of the form (1.24) 

or (1.7), depending on the kind of data m, so, under rather general conditions, 
the optimal estimator of Is also the best quadratures type formula for 
point data or for area means, as the case .may be. 

The conditions mentioned above are satisfied, for instance, when both 
ihe.geometrical covariance cov(f(P),f(Q)) and the stochastic covariance 
E ^nijOint j (P s(0t,Xj) , Q s (9it,>») are Isotropic , i.e. , functions only of 
the separation between the points p and Q . By definition of M ik the 
geometrical covariance obtained using this operator is isotropic, so C, z has 
the desired structure. A common assumption regarding good instruments 
is that the n l3 are uncorrelated, so D is diagonal. If the errors are sta¬ 
tionary, so their variances are constant, or at least constant along parallels, 
then matrix D has the required structure, and inverting C* x + D can be 
greatly expedited. In practice, however, this is not likely to be the case, 
as the number and quality of measurements will vary from region to region, 
resulting indifferent cfj both globally and along parallels. As a result, the 
best linear estimator in terms of the chosen error measure will not have the 
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quadratures form, and it will be very difficult to compute whenthe number of 
data values is very large. Nevertheless, as shown in section 3, quadrature 
formulas caa give reasonable estimates of the C a ^ with noisy data sets where 
the noise is uneven, so it would be interesting to get the best quadrature for¬ 
mula for a particular combination of signal and noise, provided that such a 
formula can be obtained without undue effort. 


2.9 The Best Quadrature Formula for non-Unlform, Uncorrelated Noise 

If the variance of the noise fluctuates along parallels, matrix D , though 
diagonal, is such that the minimization of the error measure (2.37) 

= tr [C - 2C„ F T + F(C„ + D) F T ] « * (D) (2.42) 

(see also (2.26)) 

can be very difficult with large N p , and the optimal estimator is not of the 
quadrature type. Introducing a "modified noise matrix" L , also diagonal 
and where the diagonal elements are 

(2.43) 


iuc = gj [ o?j (with k = 2N i+j) 


the following modified error measure $(L) can be defined: 

$(L) =tr [C - 2C„ F T + F(C„ + L) F T ] 

The optimal estimator for this measure is easy to obtain, and is of the 
quadratures type. 


(2.44) 


The parts of 5 (D) and * (L) that measure the sampling errors are 
identical, so any difference between the overall measures must come from 
the "noise propagation" parts tr [FDF r ] and tr [FLF ? ] . If the estimator 
(not necessarily optimal) happens to be of the quadratures type, i.e., for 
point data: 

62 - T T Xi* {sto} [f (e» ,Xj) + n t ] (2.45) 

1=0 4 = 0 

then the propagated noise is, assumming the n to be uncorrelated, 

I N-4 a N-l n N-l 8N-X 

®r?(D) =tr[FDF r ] = £ £ f £ L (X?’) 3 £ "i*j 

' dfeo tsoeo ' B = o »=° 1 = 0 4 = 0 (2. 


so 


(2.46) 


The "modified noise’,' on the other hand, is, according to (2.43) and (2.44): 





r 


fW. a N-i a*“l 

UL) =tr [FLfVI L I (Xi B ) a I (cos a m j AX + sin 2 mjAX) 

^ « = o»=° 1 = 0 j = o 


3N-1 


1 £ 7WT 

* N J = ° B=0» = 0 1=0 J=0 




(2.47) 


Comparing the expressions for [FLF T ] and for [FDF T ] , it follows 
that they are Identical, and since the "sampling" parts are also identical in 
(2.42) and (2.44), then 


tr[C -2C ei F T + F(C lt +L)F T ] = tr[C - 2C ai F T + F (C IX +D)F T ] 

(2.48) 

This means that the actual and the modified error measures must 
coincide if the estimator is of the quadratures type. 


Replacing D with L in equation (2. 38) and solving for the estimator 
matrix, one gets 

Ft * C 0 *(C IX + L)” 1 (2.49) 

where Ft. is the estimator matrix that minimizes the modified error measure 
(2.44). Because of the way L has been defined, this estimator is of the 
quadratures type, so the modified and the actual error measures coincide, as 
just shown. 

Assume that there is an estimator, different from c = f l m but also 
of the quadratures type, the estimator matrix of which is F, and such that: 

tr [C - 2C 0t F t + F(C m + D)F T ] <tr [C - 2C 0I Fj + F L (C„ + D) Fj ] 


Then, according to (2.48), 

tr [C - 2C SZ F T + F (C„ + L) F T ] <tr [C - 2C* f[ + F u (C„ + L) f[ ] 

(2. 50) 

which contradicts the fact that F L minimizes the modified error measure. 
Therefore, (2.50) cannot be trie, and F L must be the matrix of the optimal 
quadratures type estimator that minimizes the actual error measure o 2 . 

The optimal quadratures type estimator, as the name Indicates, is the 
best of a certain kind, not the absolute best. The best estimator, when no 
conditions as to its form are imposed, will not be (in general) of the quadra¬ 
tures type, unless D happens to have the "right form" specified before, 
i.e., unless D = L. 

When D = L , the optimal estimator and the best quadrature formula 
coincide. Regardless of this, the quadrature formula obtained from (2.49) 
is the best, so its error measure is a lower bound for those of all other 
quadratures formulas with the given signal and noise. 
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When D /L, while minimizing the sum of all error variances i.e., 

tr (h>), the optimal quadrature formula does not minimize each individual 
variance - To shnw fhls. ennsidnr the nrrmapalpd noise measure for 


variance a B . To show this, consider the propagated noise measure for 
when the n 4l are uncorrelated (to simplify the argument): 

0r)2» = t <x“) a r [sin} o?j = £ (xf) 2 i r o’i + if 

' 1=0 i=o^ J 1=0 J=0 i = 

The modified error measure, minimized by the formula, is 

fc-l »N-1 r 3N-1 N-l 3N-1 

I(x;*) a I&;n, i 4XJ ; I af, = l ( X ;v it Of, 

1=0 3=0 OU1 ^j-O 1=0 3 = 0 

Clearly, both are not the same, unless 


IW 


2mj AX af. = 0 


which is not lilcely to be fulfilled for arbitrary a* . However, looking at the 
reasoning which leads to (2.47), one can see that the sums of the modified and 
the actual error measures for pairs (C n „ , S nB ), and also for the individual C^ 0 , 
are alread}' identical. From this follows that the variance of the error per degree 


t t < a - <4 


a-o«=o 

and per average coefficient per degree: 


6<f„ - 


2n + 1 


(2.51) 


(2. 52) 


are also identical to the modified measure. So, while nothing can be 
predicated of individual coefficients, the error for each harmonic as a whole 
and that for the "average coefficient" in it are going to be minimum. By Parseval's 
identity (1.12), if the coefficients were used to calculate, say, geoidal undula¬ 
tions, the mean squared error of the computed geo id, globally, would be tie same 
as the sum of the error squared of the normalized coefficients, so individual 
coefficient variances are of little interest in this and similar applications,while 
the af n are very important. This shows that the optimal quadratures formula 
when D / L can be just as useful as when D = L. 

The discussion in this paragraph lias been centered on point value type 
formulas; the conclusions apply equally well to area mean type formulas, the 
extension of the reasoning being quite straightforward. 


2.10 The Structure of the Covariance Matrix and its Consequences 

The following discussion summarizes some results presented by this 
author in a previous report (Colombo, 1979a). In order to be able to calculate 
the variance of the error a® ( ® with expressions such as (2.26), and also to be 
able to obtain the optimal estimator according to collocation theory, it is neces¬ 
sary to create and invert the N p x N p matrix (C ZI + D), which can be very 
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large if the number of data N p is large. In the case of regularly sampled data 
this two problems can be greatly simplified if the covariances and the grid 
has certain symmetries. The most important of these are- (a) the sampling 

injongihidemust be at constant interval and along parallels (or parallel bands, 
i.e., rows of blocks); (b) forgiven i and p the covariances. cov(u(0i ,A,), 

(or cov^, v Pq )) and Ejn„, n p ,}, must defend 

only on ij-ql . It is also very advantageous, though not essential, that the 

grid be symmetrical with respect to the Equator. 

In what follows N r is the number of parallels and N1 the number of meridians 
(N r = N, N1 = 2N when the grid is equal angular). 


Under this set of conditions, if the data vector m is ordered according 
to (2.8) and is subdivided into partitions nij, where 

= [m, 0 m tl . . . 

includes all data values in the same parallel or row of blocks, then the 
matrix (C It +D) can be partitioned into blocks C lp , each of dimension 
N1 x N1 , containing the covariances between the data along rows i and p . 


Each block C lp has a Toeplitz circulant structure, because its elements 
satisfy the relationships 

c j« = c wq+i 5 c !o = c ri Ni“i when j > 0 


which follow from the fact that parallels are circular, and that the covariance 
between points in parallels i and p is a function of |j-q| . Moreover, 
the elements in the first row or column (the C lp are symmetrical) also 
satisfy 

c oq “ co^ x -, when q > 0 

Therefore, the first row can be represented exactly as a sum of iNl + 1 
cosines: 

= £ a‘. p cosmjjy- q (2.53) 

The a ,p form the discrete Fourier transform of the sequence 


c&. 


> »r 

■'O 2 


c 


ip 

0 Mi— 1 
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If 


r 1 , 5 = Hai ? (2.54) 

where H = |N1 if m = 0 

l Nl if m ± 0 

T 

If R(m) is the matrix where each "ip" element equals r 1 / , then (as 

explained by Colombo,(op. cit.)) Inverting (C 2I + D) is equivalent to in¬ 
verting the Nr x Nr matrices R(m) for m = 0, 1, . . . -|N1 
Isotropic covariances satisfy the "|j - qi condition" mentioned above, so, 
for fixed AX , the covariance matrix always has this regular structure. 


Let 


& - [ Q m ° 4x - {,“} m4x - • • •{r> m < Ni - i > 4x ' r 

A vector of the type 


a 


= [v 0 c 


at 

o t 


ar 

v i £.« . • 


• v Nr _ x c“ T ] 


(2.55) 


shall be called, for convenience, a vector " of frequency m" . 

Under the conditions described before ail the eigenvectors of (C„ + D) 
are vectors of frequency m , with m = 0, 1, . . . •gNl . Moreover, if 
(* = 2 » • • • Nr) Is one of the Nr eigenvalues of R(m), and if 

It. = [s 1 /. . . s«U ] T 


is the corresponding eigenvector of R(m), then X t> is also an eigenvalue of 
(C t * + D), and the pair 



is; 


t« C «T 
a » 



r 


the two corresponding eigenvectors of (C lz + D). Therefore, to decompose the 
large covariance matrix In eigenvectors and eigenvalues is equivalent to 
decomposing the \ N1 + 1 matrices R(m) , and this is why the latter are rele¬ 
vant to the inversion of (C M + D): the eigenvalues of the inverse are the 
reciprocal of the X t , , while its eigenvectors are the same as the . Further, 
this implies that (C„ +D) -1 has the same structure as the covariance matrix, 
i.e., it consists of Toeplitz-circulant blocks. 

Since (C tt + D) -1 has eigenvectors of frequency m, then, if h is 
a linear combination of vectors of a given frequency, z = (C „ + D) _1 h is 
also a linear combination of vectors of that frequency. In the case of point 
data, from expression (2.6) follows that the cross-covariances vector in 
(2.40) Is 
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T 0 = ,-TT ,, GO - r, at 

" 2nT7 11 » ( cos0 o) £= • • • P M (cos0 Mr _i) c, ] 
= [ k S“cf f . . . kg-! cf] 


(2.56) 


Define 


k 0B = [kS“. . . ksr-i]' 


(2.58) 


na r no ns 

X. L Xo • • • Xn r —i j 

* a l 

Then T n0 = (Cj, 2 + D) c_,, a, z must be of the form 

toe _ , nm a T rj! a T ., T . __ 

£ n n l Xo £ r ... Xn r-i £ a ] (2. o7) 

where, according to Colombo (ibid). 

X M - lt(m)" 1 k"” (2.58) 

Similarly, for area means, 

3 re^.e_ rx ^aa t 

£ bb0> i = -— [• . • / P D „ (cos6) sm0d9/ c“dX. . .] 

( 2n + 1 ) A *i Jdi J\ t 

= ^ . r‘* i p.„(co S e,sjnedDf! A(ra) )c,° J n(m >) ci)’. t2 ' a!>) 

(2n + 1) A jj J 01 001 VNj \ A(m)f -y 

(where A is supposed to be independent of j > according to expressions (1. 7) 
and (2. 7), so 

r fA(m)} o , jB(m) | , \* T 

ln a !• • • Xt ^ £s + | A(m| - B J * * * 1 

In conclusion, the optimal estimator for point values has the form 


Xj+AA T 


c“ dX. . .]’ 


SOL _ ,r a J .. N r 1S C _ n» fCOS\ 

£»a (£sa) £i L. L, Xi ]Sinj m jAA!Ujj 

i=0 3 =0 l / 


(2. 60) 


while that for area means is of the type 


42 - <tVa -Y Yxr I cosmJAA ♦ 


J B | m) \ sin m j AX) m 1>( 
\ (2.61) 


so they are both of th e Quadratures kind, as anticipated in the proceeding para¬ 
graph. 


2.11 Setting up and Inverting the Covariance Matrix 

Each block C 1r of (C„ + D> is wholly determined by the 1st ? N1 + 1 
elements in its fir?t row; if the number of operations required to compute any 

- 42 - 






r 


element of is k , then only $N1 + 1) k operations are needed per block, 
instead of Nl\ , as would be the case if C 1 * did not have the Toeplitz structure 
described previously. This is a reduction of the number of operations by a 
factor of 2N1 , and clearly applies not only to C 1 ® but to the whole covariance 
matrix as well. So (C tl + D) can be set up about 2N1 times faster than an 
ordinary matrix of the same size. 


The total number of elements to be computed is -|N1 x Nr 3 ", or N 3 in the 
case of equal angular grids. If the grid is a fine one, this can still be a very 
large number of covariances. This is particularly serious in the case of area 
means, because the area mean covariances are given by expressions of Ihe form 


cov(U{ j, u^,,) = M 

cov(u(0,A), u(8',\')) dodo” 


| fudo /udo , = /* r M fu(0,A) u(0',\')j 
l^lj Jo i, Jo* 1 J 


i Jo* 


dodo' 
(2. 62) 


involving double area integrals of the covariance function. Numerical quad¬ 
ratures methods, such as the one described in paragraph (4.3), have been 
used in the past to obtain cov(Ut., u^,,) (see, for example, Rapp (1977)). 
These methods take so much time in the case of fine equal angular grids 
for instance, that it may be practically impossible to use them to set up the 
covariance matrix of a global data set, in spite of the reduction by 2N1 in 
the number of operations. Fortunately, the coefficients a* p in the Fourier 
expansion of the elements 


c 


ip _ 
Ji 


cov(u 1Jf u M ) + E 



(expression (2.53)) can be obtained by means of a series expansion (truncated 
to a conveniently high degree N %x ) according to expression (4.14) in paragraph 
(4.1). These coefficients are 


where 



I i 


K Mnax 

Ia,8Nh+» tf + £ Y. In.aNh-»,l 
h=i s = a 


In»3Nh—«,p F(m) 


1 

F(m) 



Oft-A 0 

, < oos9 > sm9d6 vfrh 

AX 3 if m = 0 
(2 m" 3 ) (1-cosmAA) 


_1_ 

AX(cos 8j - cos(8j+A0)) 

and (2h + 1) N s N, tJ 


A similar reasoning to that for area means leads to an analogous formula for 
point values: 

K Na*X K NwH 

a « = [ I L (^*»3Nh+ ■ (cos 0J ) P a 3Mh»B (COS ^)) +11 (P B a N h-« (cosSj) . 
a=0 n=« - h = 0 n = « ( 2 . 64 ) 
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The importance of ( A 63) and (2. 64) is that, if the signal and noise are such 
that the number of terms in the summations is not too large (N„, * is a 
"manageable" number!, they allow the direct determination of the elements 
of the R(m) matrices according to (2.54). In this way, the R(m) can be 
created without first having to set up the whole covariance matrix and then to 
obtain the discrete Fourier transform of the first row of each C i; ' . This 
advantage further increases in the case when the grid is symmetrical with 
respect to the equator, a situation that applies to all equal angular grids. 

Then each R(m> is persyrametrical, i.e., symmetrical with respect both 
to the main diagonal and the main antidiagonal, provided D is also persvm- 
metrical (for instance, uniform noise). This means that only approximately 
-|-N r 2 elements in each R(m) are different and have to be calculated indb* 
vidually. 

Having set up the R(m) without first creating the covariance matrix, 
the inverse of (C* * + D) can be found by the equivalent operation of obtaining 
all BJmf 1 . The number of operations in a matrix inversion is usually 
0(dimension 3 ), or 0(N e, i for a covariance matrix of an equal angular data 
set. The number of operations per R(m) is 0(N?), or 0(N 3 ) for the equal 
angular grid. In fact, as explained in (Colombo, 1979a), the inversion of a 
persymmetrical R(m) is equivalent to that of two matrices of half its dimen¬ 
sion, one related to vectors of frequency m of the cosine type, and the other 
to vectors of the same frequency of the sine type. This further reduces calcu¬ 
lation by a factor of . With 0(N) R(m) matrices to be inverted, 
the total comes to 0(N 4 ) operations, or 0(N 3 ) times less than for the in ¬ 
version of (C„ + D) by ordinary techniques (Choleskii factorization, Gauss- 
Jordan elimination, etc.). 0(N rf ) is also the order of the number of data 
points in the grid, so in the case of a 1° x 1° equal angular grid with 64800 
elements the reduction in computing time is 0(64800). 

The numerical examples in section 3 all involve 5 C x 5° data sets with 
2592 elements, so (C xz + D) Is of dimension 2592. Settingup and inverting 
such a matrix is a large exercise, even with a modern digital computer such 
as the AMDHAL 470 at Ohio State, unless the matrix has a strong structure 
that can be exploited to simplify the work. As such is indeed the case here, 
the subroutine NORMAL described in Appendix B has been able to do the "’hole 
setting up and inversion in only 20 seconds. 

The inversion of (C-+ D) requires (^dimension 3 ) operations (0(N 4 )) 
instead of 0(dlmension 3 ) because of che Toeplitz clreulant structure of the 
C' 9 blocks. This "©(dimension 3 )" property is common to other algorithms 
for inverting Toeplitz-type matrices, such as the famous Trench algorithm 
(Trench, 1965), and the Justice algorithm (Justice, 1977), the first for data 
sampled on the real line and the second for data sampled on the plane, So, 
in spite of its "rather wicked" nature, the sphere allows this very convenient 
property of regular grids to apply also on its surface. In fact, not only on the 
sphere, but also on any aody of revolution (cone, oblate and prolate spheroid, 
hyjterboloids and parabotoids cf revolution, etc.,) rerJar sampling ana 
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covariances that satisfy the ” |j-q| condition" will result in covariance matrices 
of the type described here, and this is also true of other matrices based on 
symmetrical kernels, such as the normal matrices of point mass models, 
w^en the points belong to a regular grid, etc. Finally, the optimal estimator 
C nB = 5} for this type of covariance matrix is, as shown in the previous 

paragraph, of the quadratures type, so the optimal t ^ can be obtained using 
the same efficient algorithms described in paragraphs (1.5)-(1.7). Altogether, 
the powerfiil structure of the covariance matrix for regular global data sets 
is most remarkable. One of its many advantageous features is that, because 
the creation and invertion of each R(m) can be done quite independently from 
those of the others, the algorithms developed for this type of matrices are 
eminently suited for implementation in parallel processing computers. 

The separation of the algorithm according to orders also means that, 
although setting up and inverting all the R(m) may require a large number of 
o()erations, only a fraction of those actually correspond to the recovery of the 
of any given m, so the numerical errors due to rounding or truncation 
are not likely to accumulate to any great extent in the results. 


2.12 Optimal Formulas for non-Uniform, Correlated Noise 

Irregular noise, already discussed in paragraph (2.9), may be due not 
only to the varying quality of the measurements, but also to the way the data 
is "grided", L e., the way the value attributed to a node (or block) ij is 
obtained by interpolation from actual measurements nearby, as usually data is 
not sampled regularly on a global basis. As the number, disposition, and 
quality of the measurements used will vary from point to point in the grid, so 
will the accuracy of the interpolated values. Furthermore, even if the measure¬ 
ments themselves are not correlated, the grided values may be correlated 
because some of the data may be used for more than one interpolated value. 

This brings about the question of what can be done when D is neither diagonal, 
nor are the D lr blocks in D , corresponding to the Cblocks in C,* , all 
Toeplitz circulant. The answer is a simple extension cf the results already 
obtained for the uncorrelated case. 

When the noise is both non-stationary and correlated, replacing the co- 
variances E {n 14 nr,} with 

= 55 (X (E ^ n ‘ 3 nr * w ‘ ) + E {” 1J “***}>] • where h = J - s * 

will result In a modified "noise matrix" L where the (corresponding to 
the D and the C*) will be all Toeplitz circulant, because the "covariance" 

(T^^ satisfies the condition that, for a given i and r , it is a function of |j-s| 
alone. The optimal estimator for the modified measure 

$(L) = tr(C - 2C« t F T + F(Ct, + L) f} 
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The importance of (2. 63) and (2. 64) is that, if the signal and noise are such 
that the number of terms in the summations is not too large (N s , * is a 
"manageable” number}, they allow the direct determination of the elements 
of the R(m) matrices according to (2.54). In this way, the R(m) can be 
created without first having to set up the whole covariance matrix and then to 
obtain the discrete Fourier transform of the first row of each C 1 *' . This 
advantage further increases in the case when the grid is symmetrical with 
respect to the equator, a situation that applies to all equal angular grids. 

Then each R{m) is persyrametrical, i.e., symmetrical with respect both 
to the main diagonal and the main antidiagonal, provided D is also persym- 
metrical (for instance, uniform noise). This means that only approximately 
-^-N r z elements in each R(m) are different and have to be calculated indi- 
vidually. 

Having set up the R(m) without first creating the covariance matrix, 
the inverse of (C zz + D) can be found by the equivalent operation of obtaining 
all R/m)- 1 . The number of operations in a matrix inversion is usually 
0(dimension 3 ), or 0(N e } for a covariance matrix of an equal angular data 
set. The number of operations per R(m) is 0(N?), or 0(N 3 ) for the equal 
angular grid. In fact, as explained in (Colombo, 1979a), the Inversion of a 
persymmetrical R(m) is equivalent to that of two matrices of half its dimen¬ 
sion, one related to vectors of frequency m of the cosine type, and the other 
to vectors of the same frequency of the sine type. This further reduces calcu¬ 
lation by a factor of . With O(N) R(m) matrices to be Inverted, 
the total comes to 0(N*) operations, or 0(N 3 ) times less than for the in¬ 
version of (C zz + D) by ordinary techniques (Choleskii factorization, Gauss- 
Jordan elimination, etc.). 0(N a ) is also the order of the number of data 
points in the grid, so in the case of a 1° x 1° equal angular grid with 64SOO 
elements the reduction in computing time is 0(64800). 

The numerical examples In section 3 all involve 5° x 5° data sets with 
2592 elements, so (C zz + D) is of dimension 2592. Settingup and inverting 
such a matrix Is a large exercise, even with a modem digital computer such 
as the AMDHAL 470 at Ohio State, unless the matrix has a strong structure 
that can be exploited to simplify the work. As such is indeed the case here, 
the subroutine NORMAL described in Appendix B has been able to do the "hole 
setting up and Inversion in only 20 seconds. 

The inversion of (C- z + D) requires (^dimension 3 ) operations (0(N 4 )) 
instead of 0(dimer.sion 3 ) because of ihe Toeplitz clrcuiant structure of the 
C 15 blocks. This "0(dimension a )" property is common to other algorithms 
for inverting Toeplltz-type matrices, such as the famous Trench algorithm 
(Trench, 1965), and the Justice algorithm (Justice, 1977), the first for data 
sampled on the real line and the second for data sampled on the plane, So, 
in spite of its "rather wicked" nature, the sphere allows this very convenient 
property of regular grids to apply also on its surface. In fact, not only on the 
sphere, but also on any oody of revolution (cone, oblate and prolate spheroid, 
hyjierboloids and paraboloids of revolution, etc.,) rer.lar sampling and 
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covariances that satisfy the " j j—q| condition" will result in covariance matrices 
of the type described here, and this is also true of other matrices based on 
symmetrical kernels, such as the normal matrices of point mass models, 
wj^en the points belong to a regular grid, etc. Finally, the optimal estimator 
C a> = J^am ® for this type of covariance matrix is, as shown in the previous 
paragraph, of the quadratures type, so the optimal 6^ can be obtained using 
the same efficient algorithms described in paragraphs (1.5)-(1.7). Altogether, 
the powerful structure of the covariance matrix for regular global data sets 
is most remarkable. One of its many advantageous features is that, because 
the creation and invertion of each R(m) can be done quite independently from 
those of the others, the algorithms developed for this type of matrices are 
eminently suited for implementation in parallel processing computers. 

The separation of the algorithm according to orders also means that, 
although setting up and inverting all the R(m) may require a large number of 
operations, only a fraction of those actually correspond to the recovery of the 
Cn,, of any given m, so the numerical errors due to rounding or truncation 
are not likely to accumulate to any great extent in the results. 


2.12 Optimal Formulas for non-Uniform, Correlated Noise 

Irregular noise, already discussed in paragraph (2.9), may be due not 
only to the varying quality of the measurements, but also to the way the data 
is "grided", L e., the way the value attributed to a node (or block) ij is 
obtained by interpolation from actual measurements nearby, as usually data is 
not sampled regularly on a global basis. As the number, disposition, and 
quality of the measurements used will vary from point to point in the grid, so 
will the accuracy of the interpolated values. Furthermore, even if the measure¬ 
ments themselves are not correlated, the grided values may be correlated 
because some of the data may be used for more than one interpolated value. 

This brings about the question of what can be done when D is neither diagonal, 
nor are the D lr blocks in D , corresponding to the C lr blocks in C lz , all 
Toeplitz circulant. The answer is a simple extension cf the results already 
obtained for the uncorrelated case. 

When the noise is both non-stationary and correlated, replacing the co- 
variances E {njj n,,} with 

5>rW = 4S (% (E f n ” *** ) + E { n n °ri-» } >) * where h " J ‘ s » 

will result in a modified "noise matrix" L where the i/ (corresponding to 
the and the C*) will be all Toeplitz circulant, because the "covariance" 
CT Sr|h| satisfies the condition that, for a given i and r , it is a function of |j-s| 
alone. The optimal estimator for the modified measure 

«(L)=tr(C - 2Cat F T + F(C z , + L) F^ 


- 46 - 



must be of the quadratures type, because of the structure of L . To show 
that is is also the best estimator of this kind in terms of the original norm 


#(D)= trf C - 2C ei F T + F(C It + D) F*| , 
the proof will proceed much as in the case of paragraph (2.9). 

The propagated error measure for C 3B is 

( N-i ai+-i n-i a t*-i , •. i 

<&. = E {| 0 M mi 4 X ”",l 0 .L x '”{ r 4 rasiXn 4 


M_1 3N-1 N-l a«*-i (cos) (cos) 

-I £ £ £ * W m ^ x W 


msAX E 


1=0 j=Or=o*=0 
so, for C n , and combined. 


(n u n r ,J 




N-l a N-l N-l 3N-1 


I CT rr“ 

or=o 


£ I £ £ cos m(j-s) AX E |n 13 n r .} 

1=0 3 = 0 r = 0 •= 0 


thus 


N-i, n N-i N-l «N-l aN-i 

ijjlD) = tr [FDF T ] = n £ £ £ X Xi‘Xr“ £ £ cosm(j-s)AX E{n n iv4 

M_1 a N-l N-l — 0 1 ~ N-l 3 _ Sn-1 0 

= £ £ £ £xi*Xr“ £ cosmhAX £ (E{n l3 + E \% j n, jj,J) 

a=0i=01=or=0 h = ° 3 = 0 

Replacing both E [n 13 n,j*,} and E in the last expression with 

a I^,ll, is the same as replacing D with L , so 




T a N-l M-l N-l 3N-J 

4 (L) = tr[FLF]=£ £ £ £ Xi" Xr £ cosmhAX £2a 1 

v a = 0 a —o 1 = 0 r%0 l4o 1 = 0 

"r* f ^ ^ a. . .. _.._ 1 r | h | 

= £ £ £ £ Xi XT £ cosmhAX 4No 

a-0 « = o i — o a-o 

N-l a N-l N-l N-l 3 N-l . . . 

= £ £ £ £ Xi Xr £ cosmhAX £ (E{&u d,,* ) + Eln^n,^}) 

a = 0 »=0 H ° r =0 3-0 

because of the definition of ? Ir iil . Comparing the expressions for ^(L) 
and for ? (D) it is clear that they are identical. From this follows that the 
modified error measure $(L) coincides with $(D) when the estimator is the 
optimal estimator of the quadratures type for $(L), and that this must be tiie 

optimal estimator of the quadratures type for $(D) as well. The other con¬ 

clusions arrived at in paragraph (2.9) for the uncorrelated case apply equally 
well here. 


2.13 Least Squares Adjustment, and Least Squares Collocation 
(a) Band-Limited Signal 

If there is a degree N lt , above which the degree variances a B are all 
negligible or zero, then the signal can be said to be band limited, and the data 
will satisfy equations of the type 
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(2.65) 


N#»x n I _ n _q; 

m 13 = r f Y c“ y an (6,.A,) + n ti 

n = 0 or: O q£zz 0 

(tlie treatment here is for point values; the extension to area means is trivial) 
With one equation such as (2.G5) per point in the grid, the result is a system of 
equations 

m + v = Ac (v = -n) (2.66) 

where A is a N p xN, matrix (N p is the number of points in the grid, 

N 0 the number of coefficients). The columns of A consist of elements of the 

type 


a" = y,.(e t .x J ) 


(2. 67) 


According to the discussion in paragraph (1.3), if the grid is equal 
angular, A has full rank when N,. x * N . N c = N 2 and 
the upper limit in the summations is N - 1 in what follows. 


Least squares adjustment is a method for solving for the while mini¬ 
mizing the propaga te d noise defined in paragraph (2.4). The least .squares 
solution is 


c = (A T D -1 A) _1 a’d" 1 m 
= G _1 A T D' 1 m 


( 2 . 68 ) 


where 

G = A t D -1 A (2.69) 

is the Nc x Nc normal matrix, while D = E{n n 1 } is the same noise matrix 
considered before. Clearly, the least squares estimator matrix is 

= (A t D _1 A) A T D _1 

When the noise has zero moan (E {n} = o ), the estimator of (2.68)is the best 
linear unbiased estimator, because it minimizes 

t r (E {Fn.n T i,,T } 1 = tr [FDF t ] 

and E {F {* >- n)} = £ • If, in addition to all this, the probability distribution 
of the noise is Gaussian, then (2. 68) corresponds to the maximum likelihood 
estimator as well. In many scientific applications the noise has approxinntely 
zero mean and near-Gaussian distribution, while D is known reasonably well; 
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for this reason, methods based on expression (2.68) are used quite often. The 
linearity of the resulting estimators is helpful, because this avoids the use of 
methods based on non-linear formulas that are usually difficult both from a 
theoretical and from a practical point of view. The variances of the estimates 
are given by the corresponding diagonal elements of the a posteriori variance- 
covariance matrix 

E = (A T D" 1 Af 1 = G" 1 

It 

Therefore, to obtain both the estimates and their variances it is necessary to 
know CT 1 . Sometimes, because of the nature of A and D , G can be seriously 
ill-conditioned, the Inversion suffering from strong numerical instabilities. To 
reduce this problem,a simple device known as regularization is often used (see, 
for instance, Tikhonov and Arsenin, 1977). Generally speaking, regularization 
is the introduction of a slight change in a problem, so the solution virtually 
rrmains the same, but the modified problem has better numerical properties. 

In least squares methods regularization usually implies adding a small positive 
definite matrix K(diagonal, as a rule) to G before attempting to invert it. 

The regularized optimal estimator would be 

c = (A r D' x A + K) A T D“ l m (2.70) 

The inverse of the covariance matrix of the harmonic coefficients C is a pos¬ 
itive, diagonal matrix which could be used to regularize the normal matrix: 

4 = (A t D _ 1 A + C' 1 )A T D' 1 m (2.71) 

It is easy to see that this expression minimizes the quadratic form 

Q = c T C -1 c + v T D -1 v (2.72) 

subject to the constraint 

m = Ac + v (2.73) 

Moreover, (2.72) is the equivalent of the least squares collocation error measure 
when then signal is band-limited (Moritz, (1980)). This idea has been used, among 
others, by Schwarz (1975) for the determination of low degree zonal coefficients 
of the geopotential, and by Lerch et al.(1979), who employed it to estabiiize the 
adjustment of the GEM-9 gravity field model with remarkable success. The 
equivalence of (2.71) to the collocation estimator is true only for band-limited 
signals; to the "real world" the gravity field has infinite bandwith, so (2.71) 
is no more than an approximation. The band-limited assumption is a reasonable 
one, however, as the eventually become negligible for large n.. This is 
particularly true at satellite altitudes; in any case, geodesy is a science of 
wise approximations. Moritz (1980) has provided a very clear and concise 
explanation of the use of collocation in general, and expression (2.71) in parti¬ 
cular, in spherical harmonic analysis. 
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An alternative derivation of (2.71) follows from the matrix equation 

CA t (ACA T + Df 1 - (A t D” 1 A + C" 1 )' 1 A t D _1 (2.74) 

(see, for instance, Uotila (1978), equation (29)), which is valid for symmet¬ 
rical matrices, provided the inverses or pseudo inverses of D , C , and 
(A C A t + D) do exist. According to the definitions in paragraph (2.2): 

C 0 * - m{cz’} = m {££. TaT } = m {££ T } aT = caT ( 2 * 75 ) 

C„ = m{z z T } = m{ac c T A T } = Am{c c T } A T = ACA T (2.76) 

Replacing C 0 * and C tt in the expression of the collocation estimator matrix 
(2,39) with their equivalents given by (2.75) and (2.76): 

F = C„ (C„ + Df 1 = CA T (ACA T + Df 1 (2 .77) 

= (A t d" 1 a + c " 1 )' 1 a t d _1 

according to equation (2.74). This shows that the "regularized" estimator matrix 
(A t D -1 A + C -1 ) -1 A t D -1 Is indeed the same as the collocation estimator 
matrix, so(Z, 71) represents an alternative form of collocation when the data 
is band-limited. 

(b) Infinite Bandwith 


In this case the "observation equations" are 



00 

m 13 = £ 

o = 0 

a 1 /v 

I E 

■= 0 a-o 

Yf,(8i ,X 3 ) + n„ 


Calling 

» a 1 « 

- l E E 

a=H*~L a = 0 

*?.<«t.x,) 

(2.78a) 

and 

W - [w x . 

. . w^ . . . w Np ] r (k = 2N 1+ j) 

(2.78b) 

then 

N 

m u “ I 

a J. -a 

E E c » 

( 0 i *Xj) + w, 3 + n u 

(2.79) 


a =0 iso 0 


and, regarding this expression as a modified observation equation, and 
replacing D with D + m|ww t | in (2.71), the linear estimator that mini¬ 
mizes the quadratic form 



(2.80) 


£ = (A T (D + m{w w t } f 1 A + C -1 l" 1 A T (D + M|w w T } J" 1 m 

It is easy to show, either following the lines of Moritz (1972), or going 

-50- 


(2.81) 




back once more to the matrix identify. 74), that expression (2.81) is identical 
with the estimator of least squares collocation. Rigorously speaking, expres¬ 
sion (2.81) should be used whenever of ^ 0 for n > N , even if = 0 for 
n > N,, x t for some finite N # x > N . 


2.14 Ridge Regression and Least Squares Collocation 
Consider once more the estimator 
c - Fjjs m = (A T D _1 A) -1 A T D~ 1 m 
If there is no noise, so m = z_ = A £ , and if n < N , then 
c * (A T D -1 A) -1 A t D _1 A c = c 

According to the definition given in paragraph (2.4), the sampling error of 
this estimator is zero, so the measure of this error must be also zero. If. 
the noise has zero mean, it follows that 

E{6} =■ E{ Fjt . m} - F ( . E {a c + »} - F,. A c + E {n} = c 

or, as it is usually said, the estimator is unbiased. Moreover, by a simple 
extension of the Gauss-Markov theorem to the case of a general symmetrical 
positive matrix D (see, for instance, Bibby and Toutenberg, 1977), F^, m, 

of all linear vmbiasecLestlmators, has the least propagated error measure 
trjFjj, DF T £„j = tr - (A T D” l A)j , as mentioned previously. 

The estimator of expression (2.77) does not, in general, give perfect 
estimates of c in the absence of noise: it is a biased estimator, and the 
measure of the bias is tr |c - 2C„ t F T + FC IZ F T j (thisternf’ean no longer 
be regarded as the measure of the sampling error, as it is the presence of 
CT 1 inside the parenthesis in (2.71) and not the sampling that brings about 
this error). According to (2.39), F is the estimator that minimizes the 
total error measure, so 

of = tr {c-2C sj F 7 + F(C„ + D) F T } 

< tr {c - 2C, 2 F T x , + F*. (C„ + D) F^,} 

- tr fo.Dl',.} 

If the covariance matrix is positive definite (which only requires that all 
/ 0 for 0 < n r N) then the last expression applies with'strict inequality 

tr {c - 2C., F T + F(C,, + D) F T } <tr{F/, DfJ.} 

- tr{ (A T D _1 A)" 1 } 
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In the band-limited case:. 





Some may find this result rather surprising: the best estimator with zero bias 
is in fact worse than the biased estimator of (2. 71) ! The difficulty is only 
apparent: F^. is the best estimator with no bias; once the "no bias" con¬ 
dition is removed, the expression above merely indicates that there is an 
estimator in the larger class of the estimators that have a bias (including 
those with zero bias) such that the sum of the bias and the propagated noise 
is smaller. From this, it is clear that 

tr |FDF T | < tr j(A T D“ 1 A )| (as tr Ic - 2C« F T + FC„ F T j 2 0 ) 

which is indeed possible when the condition C - 2 C 0l F T + F C« F T =0 
is removed. In fact, there is nothing very new about all this: the use of 
biased estimators to obtain estimates with small variances is a reasonably 
well-established practice in applied statistics. In particular, the technique 
known as ridge regression consists in using the biased estimator 

c = (xx + K) _1 x T m 

with a suitable choice of K (Bibby, 1972). Clearly, this expression is the 
same as (2. 71) when x = A, D=I, and K = C -1 . With some obvious modifi¬ 
cations suggested by (2.81), this argument can be extended to the estima¬ 
tion of c when n a N , so it can be said that within the scope of spherical 
harmonic analysis least squares collocation is a form of ridge regression . 

This brings up the question of just how realistic the error measure 
is; after all the best of all possible estimators in terms of a given norm 
could be a very had one for some specific problem where that norm is not 
suitable. To answer this question, one must start by defining the meaning of 
"realistic". If one is interested in minimizing the actual error variance of 
the coefficients per degree, i.e., the expression 

6 » = E t A - Cf. f (2n + If 1 
a= o i=o 

which may be of interest because this corresponds to, say, the global n mefm 
square of the error of representing thp continuous function (9,A)= £ 

^nn T n ,(0,A) with f (B , X)~ q^ o ^»» (® »^), according ~ 
to Parseval's identity (exp. (1.12)), then one could say that a realistic 
measure is one that gives close estimates of the actual crrorvariances. 

The "actual error" measure 6® corresponds to one of infinitely many 
’fevents" over which the collocation measure is an average. The proof of a 
pudding being in the eating, the reader can judge just how realistic the collo¬ 
cation measure is by looking at the numerical results in section 3, where the 





value of the measure turns out to differ only by a small percentage from the 
actual variance in each one of a number of simulated "events", j.e., the 
recovery of the Cj^ from"simulated data',' where the q® are known random num 
bers scaled to have the desired power spectrum. 


2.15 Structure of the Normal Matrix 


The elements of the normal matrix G (in the case of point data) are 
of the form 

N-l 3N-1 

Sms’,~._X X 

N-l aN-i ^ . (2.82] 

= X P nB (cos0j) Pp^cos0j) X \s?n} m J AX ■isi > n) qjAXor "i a 

where for all 0 <j ^ 2N (i.e., "regular" noise as defined in 

paragraph 2.8). If n , m < N , then the following equations apply: 

Tte} m J sx {S}'!i ix = o «< 2 ' 6s 

3 =o m f q 

Moreover, if the grid is symmetrical with respect to the Equator (as equal 
angular grids are), and if of = af-j-i (for instance, if is indepen¬ 
dent also of i ) the relationships 

of P n , (cos©,) = P nB (cose;) (-l)"""ojj5- 1 

( cos0 i)=P Pa (cos Q\) (-l) ! ’~ q ctn^L { 
ei = 77-e t 

must apply, according to par. (1. 2), and from these follows 

N_1 _ 

|T erf 3 Fjj^osSt) P Pb (cos8j) =0 if n - p is odd. (2.84 

In brief: if n , m < N , and the grid is regularly spaced in longitude, and 
symmetrical with respect to the equator, then 

la f 0, m 4 q 


« 0 if 


(2.85 


n - p is odd 


If neither of the conditions listed above apply, then g BB( p, may or may not 
be aero. If the coefficients C^, are ordered in c so that all those of the 
same m are grouped together, and for a given m all "C BB are separated 
from all C* B * and further more all C®, with n - m even are separated from 
all those with n - m odd, then the normal matrix becomes arranged in such a 
way that all potentially non-zero elements (i.e., not satisfying (2. 85)) are also 

A 

grouped together forming a series of diagonal blocks G„ , where 6 signals 





the parity of n-p. Each one of these diagonal blocks is made 
exclusively of one of the following types of elements 


»"P even (6=0) 
gX » n-p odd (6= l) 


If the .grid is not symmetrical with respect to the equator, groups of type 

Jand \§=i J , or IqJ and become included into larger non¬ 

zero blocks, as there are more non-zero elements in that case. However, 
here the discussion will cover only the symmetrical case. 


The largest blocks are Go'^, and their dimension is N 2 ; 
the smallest blocks have dimension lxl, for example . 


The inverse of any block-diagonal matrix suoh as G is another block- 
diagonal matrix made up of the inverse of the blocks of G . There are 
4 N - 2 diagonal blocks G^’ in G , and as many in G* 1 . 


The eigenvalues of G are those of the diagonal blocks, and the 
eigenvectors of G , those of the same blocks ’'expanded" with zeroes at both 
ends, so as to reach the dimension N p of G . 


The estimates' vector 


c = G -1 A T D" 1 m 


A 6 

can be partitioned in the same way as £ , each partition c®’ including all 
coefficients' estimates of the same m , oc, and parity 6 of n - m ; 


*a, 6 
£a 





D -1 m 


( 2 . 86 ) 


where A B isa(N - m) x N p matrix with rows that are N p - vectors 
of the type 


<*«6 , - « rcos-< 

a»,n = t • • • Pn.fcosej) | sin | mjAX. . . ) (2.87) 

(n - m even if 6=0, odd if 6 = 1 ) So the rows of (G^’ ) _1 are 

linear combinations of vectors of the same frequency m and, therefore, also 
vectors of the same frequency: 


ILmi = [ • . • Xi" jfehf} mj AA. . . ] (2.88) 

Consequently, the estimate of a giver C nB is 

6“ = h“ D^m = *? Z Xi“ {shi } m 3 AA of mjj (2.89) 

i = o 3 = 0 

and this is a quadratures type estimator. 
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Setting-up and inverting the G, blocks is tantamount to setting up 
and inverting the whole matrix G . Since all operations related to one of the 
blocks are independent from those for the others, the inversion of G i|^ g 
ideally suited for parallel-processing computing. On average, each G a ' 
requires 0(N 3 ) operations to invert, or 0(N 4 ) altogether. Inverting G 
by ordinary techniques would involve 0(/(N ) 3 ) 3 ) = 0(N 6 ) , so there is an 
increase in efficiency of 0(N 2 ) . 


Finally, it is quite simple to show that these properties carry over both 
to the case of area means, and to problems where the surface being studied is 
not a sphere, but a surface of revolution symmetrical about a plane perpendi¬ 
cular to the axis of rotation, provided that the longitude increments be constant 
and the grid symmetrical with respect to the "equator". In this lattercase, 
the expansion of the signal in solid spherical harmonics is 


i] 


N-l 


an~\ 

C 

isO 


E 

a=o 


a 

“ 5 TT 


— a 
Y 


(®i > Xj) 


and the factors r 1Ui are symmetrical about the equator, from which all the 
properties already mentioned for G follow. 


Clearly, the structure of G possesses many properties similar or 
identical to those of (C lz + D) when the data is regularly sampled on a sur¬ 
face of revolution. These similitudes underline the intimate relationship 
between least squares adjustment and least squares collocation shown in the 
preceeding paragraph, to toct, as least squares, regularized least squares, 
and collocation differ only in the diagonal matrix ( K , orC" 1 ) being added 
to G in expressions (2.68), (2.70), and (2. 77), all the properties mentioned 
here for G apply to the normal matrices in each of the throe methods equally 
well. The one important consideration, in the case of collocation, is that tj»e 
data be band-limited . Otherwise, expression (2.81) indicates that (C+M’;Wwm 
and not C* 1 mu3t be added to G . Matrix C + M-i w w T \ has the same 
Toeplltz-type structure of (C zz + D) discussed in paragraph (2:! 10). Therefore, 
creating and inverting the normal matrix requires: (a) creating and inverting 
C + M-fw w r |, and (b) creating and inverting A T D -1 A + (M jw v/'l- + c , 
which can be shown to have the same block diagonal structure discussed here. 
This is twice the work needed to set up and Invert (C z z + D) using the approach 
of paragraph (2.11), so, in the case of infinite bandwith, that approach is more 
economical in computing and, therefore, more practical. 


There may be one important point in favor of using formula (2. 77) 
or (2. 81) rather than formula (2.39) for obtaining the optimal estimator 
matrix F at least in the band-limited case: as the density of the grid in¬ 
creases, matrix (C* + D) becomes increasingly more ill-conditioned, because 
the closer distance between data points results in covariances that have much 
the same values in consecutive rows or columns. On the other hand, the non-zero 
diagonal blocks in G are likely to become more and more diagonal-dominant 
as A0, AX - 0 . This will depend on D : for instance, if the variances of the 
noise were of the form 
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then 


ofj = sin 8 t 


Li, . 

A 9 AX * a _ 6 _ 

N-l 

= Ll £ 

A0, AX 

4 tt 

A9, AX 1=0 

-0 -*0 

= -1_ 

f Y?.(e,X) 

-0 -0 

(e,X) do = 

4 tt 

Jo 


T&Oi.XO Yp^St.Xj) sine, A6AX 
3 = o 4 tt 

0 if n t p 
1 if n = p 


because of the orthogonality relationships. In general, the variances of the 
noise are not going to follow a sinusoidal law, but one may reasonably expect 
(at least with more or less homogeneous noise) that the stability of the normal 
equations will not deteriorate with A 6, A X - 0 . 


2.16 Global Adjustment and Collocation with Scattered Data 

The efficient set up and inversion of the covariance matrix (C„ + D), 
or of the normal matrix G , depend on the regular nature of the grid. If 
not all nodes or blocks in the grid have data associated with them, the data is 
said to be scattered. The blanks or "holes" in the grid destroy the orderly 
structure of the matrices, making the application of the techniques previously 
discussed impossible. Yet so strong is this structure that, even in fragments, 
still it can be dealt with more efficiently than in the case of ordinary matrices 
of the same size. 

(a) Full Region Bound by Lines of Latitude and Longitude 

In the case when there is data at every point or block inside a "square" 
region limited by parallels and meridians, the partitioning of the data vector 
along the arcs or parallels inside the zone reveals a strong structure in the 
(C,, + D) matrix, if all the other assumptions made in paragraph (2.10) still 
apply. 

If N r is the number of rows and N, the number of meridians that 
cross the region, then the covariance matrix will consist of N? blocks C lp of 
dimension N« , both per symmetrical and Toeplitz, though not circulant 
(i.e., the relationship = cjj lq+1 is fulfilled, but not c 0 * = c N ^,_i); 
moreover the first row in each block does not have the property that 
Cq, =c 0 £ 0 _, . clearly, though weaker than in the case of a global grid, 
there is a definite structure here that can be exploited to make both setting 
up and Inverting the matrix more efficient. 

Because each block C * is Toeplitz, only the N 0 elements in its 
first row have to be computed, or about £n« NT for the matrix as a 
whole, instead of ^N, a N r 3 ; this amounts to a reduction in operations by 
a factor of N, . 
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The solution of the equation ! na = (C« + D)" 1 c, for the optimal 
estimator vector f for C, ia , can be obtained by a technique such as 
conjugate gradients, or similar, in which a finite number K of matrix- 
vector multiplications (C„ + D)v t (where v 0 , v* , ... v t .... v.. are 
K intermediate N P - vectors created during the solution) constitutcTthe bulk 
of the computing effort. There is no need to go into the details of any specific- 
technique, as the reader will find excellent descriptions in the literature 
(Householder, 1964, Luenberger, I960). A discussion of the matrix-vector 
operation is sufficient here. 


Let m 1 be the N c - vector partition of the N* N.data vector m , 
containing the measurements along the ith parallel in the region, and let v/ 
be the corresponding partition in any of the v t vectors. The product (C Z2 + D)v t 
= p t is, under such partition, 

£t =[£?.. . £t T ~Y 


with 


Nr_l 




( 2 . 90 ) 


so the whole matrix-vector multiplication can be broken up into Nr products 
C Ip v^ = h lp . Because C ip is a Toeplitz matrix, the N 0 components of 
h t lp can be obtained by "weighted running averages" or discrete convolution 
of the elements of vf with those in the first row of Such convolution 

can be calculated efficiently using the Fast Fourier Transform algorithm (see, 
for instance, Brigham, 1974, Ch. 13). Therefore, all N? products involve 
0(N C N^) operations, and since there are K matrix-vector multiplications 
in the whole procedure, the total number of operations needed to obtain i ® 
amounts to 0(KN„ N, 2 ). For conjugate gradients, K does net exceed (in 
theory) N p =N r N 0 , so there should be 0(N? N 3 ) operations altogether. 

If (C z j + D) were handled by conventional techniques, disregarding its well 
defined structure, the number would be 0(N„ 3 N 3 ) , so the increase in ef¬ 
ficiency is 0(N„) , the same as for the setting up. 


(b) Arbitrarily Scattered Data 


It is common in geodesy and in geophysics to have a set of measurements 
scattered throughout the globe, without the data being on the nodes of a reg- 
ular grid or without all o"j j being equal along parallels (nonhomogenous noise). 
If the set is dense enough, however, it is possible to interpolate the data quite 
reliably on the closest nodes of a conveniently chosen grid. Assuming that this 
is done, and that the accuracies of the interpolated values are 'mown well 
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enough, then the problem can be dealt with by conventional least squares or 
by collocation. In general, there will be blanks or "holes" irregularly dis¬ 
tributed over the sphere, the actuid data points falling among them in no 
precise pattern. This problem will be considered hero as a least squares 
adjustment problem. If the data has little or no power above the Nyqulst 
frequency of the grid on which it has been interpolated, then the extension of 
the ideas that follow to collocation is quite simple, according to paragraph 
(2.13). 


Consider the element g a ^ 4 of the matrix G - A T D 1 A 


“’0 = 
bn mJ>Q 


N~1 


£ P n8 (COS0J) p p4 (cos Oj) £ 
1=0 1=0 


|sin | |sm) 

(2.91) 


wtere Wjj = 0 if the point ij is blank, otherwise YV^ = 1 ; ofj is the var¬ 
iance of the noise njj . In general (Ttj is a function both of i and of j . 
Clearly, matrix D is taken to be diagonal (i.e., uncorrelated noise). 


From the relationships 


cos m j AX cospjAX = i [cos (m + p) j A X + cos (m - p) j A X] 
cosmjAX sinpjAX = i[sin(m f-p) j A X-sin(m - p) j AX ] 
sinmjAX sinpjAX = -| [cos(m + p) i AX-cos(m - p) j A X] 


and calling 


au-i 


Cr 01 = 


c‘°= 


cos r JAXcrTl w, : 


8N-1 


Cr 1= 2 £^ sin r j A Xo,j W u 
where -N <r<2N , follows 


2N-1 

1 ^-o 


»?, W„H n*U ffipW - tc;V + 

(/3) 


(GO 

Moreover, because 


cos (-r) j A X = cos r j AX 

sin(-r)jAX= -sinrjAX 

cos r j A X = cos (2N - r) j A X (2 nA X= 2u) 

sinrjAX = -sia(2N - r) j AX 


(2.92) 


(2.93) 


follows 


o! r « - 

c, ia 

(-i f 

- 

varr “ 

c r l « 

(-i f 
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where the last one is of special Interest when r s N . So only the C ,‘ a 
with 0 < r < N are needed. Finally, calling 

= P a .(cos 0,) (cos 9j) (CL?" P| + (-D 30£ * a ] (2.94) 

(where jot- /3| is the absolute value of a- $) 


It is 


& a» # pq 


v-1 


= E aittf 


1 = 0 


(2.95) 


Assume that out of N rows the grid has only I with any data in them, so 
at least one ^ 0 in each. Once the corresponding C r lQf have been ob¬ 
tained by computing the discrete Fourier transform of o^ 3 W t j along each 
row with data (C^N 3 ) operations for the whole grid), what remains is to get 
the I non-zero terms that form each element g ffi {? a of G . As there 

are about such elements that are different, and I is O(N) (except for 

very sparse data sets), the total number of operations needed to create G 
is OfN 3 ) + 0(N & ), or virtually 0(N 6 ) . If the were computed according 

to (2.91) as it is written, instead of according to its reduced version. (2.95), 
the number of operations would be 0(2N 3 (N) 4 ) , or 0(N a ), so the gain in 
efficiency allowed by this approach Is O(N) , which is the same as in the case 
studied in the first part of this paragraph, where the data completely filled 
a "square" sector of the sphere. This count does not include the time needed 
to obtain the (cos 9 t ), as these can be pre-computed once and kept on 
disk or tape for repeated use. 


The number of operations can be reduced further, almost by half, by 
takiM advantage of the fact that F a ,(cos £ i) is a common factor in all 

w tth the same n and the same m . Furthermore, if there are pairs 
of rows i and N-l-1 (l.e., symmetrical with respect to the Equator) where 
both rows contain some data, then 

Safcpq + &na,P9 ~ ^a» (COS 9j) Ppg (COS 0 j) 

dciir'* 1 *(-1)"^“^') - *(-D" ci>“- 81 ,J 

(2.96) 


which leads to further savings in computing at f he cost of additional program¬ 
ming complexity. These economies are important, but they will not bring the 
number of operations much below 0(N S ) unless the grid is so sparse that 
I is much smaller than N . 

Notice that the normal matrix G = A T D~ x A is created here without 
actually forming the observation equations matrix A . This means consid¬ 
erable savings in computing and in storage requirements. As for the right 
hand side of the normals A T D" 1 m= b , the elements of b are given by the 
formula 
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(2.97) 


r 


N-l 3N-1 j- i 

bk = i?o 3 ^o Pa * (COS9l) IsSiJ m 3 AXor » w„ rni, 

where k=n 3 +an + m + l and W t s ■ 0 if there is no data at the point lj , 
as before. Such expression is of the quadratures type, and can be computed ef¬ 
ficiently by the corresponding algorithm of section 1, also without first creating 
A . Finally, the residuals vector v = m - A A » usually of interest, can be 
obtained as the difference between the data m and the values of 

I t (cos9 t ) mjAX 

Ct = 0a = o »“q 

computed by means of the appropriate synthesis procedure given in section 1. 
Thus v can be found without knowing A explicitly. 


The normal equations can be solved by means of conjugate gradients 
or a similar method involving M matrix-vector products x t = Gh t , where 
h t is part of a sequence of Intermediate vectors h 0 , hi . h t , . . . h, . 

Introducing the notation 

a 

fea«t = ikt 


where k=n a + an + m + l, and calling G 1 to the matrix of all g 1 ,®^, , 

SO G = £ G and x t = [ G h t = £ x t (2.98) 

1=0 1 = 0 S =0 

then each element t of the product vector is of the form 

- f 1 z t *'“£ > 6 « - p„(coss,i j- [C 'Jz< 

P= 0 q =0 0=0 Pq0 


. P p q (cos 9 t ) h^, t 


(2.99) 


Every h^ t above multiplies always the same F p , (cos9 t ); calling 
dp^ 3 P p9 (cos 9 t ) hp^ t 


expression (2.99) becomes 


*£ = Z.icoseo Z [C'M-irCl dp? 


a'-B 


P qjS 


All quantities inside the square brackets are constant if oc, /3, m , and q 
are the same. Grouping equal factors together 

.Ml ■ ' ^ 


= P». (cos 90 £* £ [cjT + (-l) 3a+8 C^l Z dpf 

<i=b!J*o p = o 

- P„ (COB 9,) £■ £ of 

where D* - £ ’ “ ** 
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1 


There are (N 3 ) products needed to form all ; (N 2 ) sums to compute the 
2N , and (N) 2 further multiplications bv the P riB (cos 0,) to form all 

As there are I non-zero G 1 , there arc, per matrix-vector product, 
IxC^N 3 ) operations, or 0(N a ) if the data is not too sparse. M , the number of 
matrix-vector products in the solution, is ©(N 3 ) , so the total comes to 0(N b ) . 
Inverting G by the usual methods and without having regard for its structure 
involves OfN 6 ) operations. The gain inefficiency is, again, O(N) . 

Besides providing a convenient way of demonstrating how the properties 
of G can be exploited to make its inversion more efficient (or the solution 
of the normal equations, both approaches are equivalent), conjugate gradients 
is interesting on its own right. Sparse data sets with a poor distribution will 
result in ill-conditioned normals, so the inversion of G may be numerically 
impossible. So-called iterative methods, such as conjugate gradients, usually 
improve the initial guess (represented by hp) of the correct values of the un¬ 
knowns, at least for the first few iterations. Improvement here means a reduc¬ 
tion in the quadratic form being minimized, such as the mean square value of 
the residuals. If the initial guess is a good one, and present day spherical 
harmonic coefficients of the gravity field are reasonable good for degrees up to 
30 or so, then a few iterations are likely to improve this guess,and to produce 
reasonable estimates of those coefficients that,being wholly unknown, are 
taken to be zero at the start. If a "few” iterations are much less than the max¬ 
imum N 3 , then a reduction in computing time of 0(N a ) takes place. This 
might allow scientists to "extend" existing models to much higher degree end 
order than at present, simply by obtaining approximate solutions of this type. 

Clearly, all that has been said here regarding G and least squares 
adjustment applies equally well to (G + C l ) and least squares collocation. 
Though the formulas have been developed on the basis of a point values' 
formulation, their extension to area means is not difficult. 
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While the descriptions of the methods for estimating the Sf, 
discussed here have been confined to the case where all the data are of one 
kind (i.e., gravity anomalies only, or magnetic anomalies only, or geopo- 
tentlal numbers • etc.) their extension to mixed data sets is immediate, 
provided that all values are globally distributed according to grids of the same 
A X, although the latitudes need not be the same as well. 

Note on the Accuracies of the when Using Conjugate Gradients ; 

Besides the estimates of the coefficients, one usually wants to know the 
accuracies of those estimates. When a few iterations of conjugate gradients 
are used to update some initial estimates in the efficient way described above, 
what are the accuracies of the improved coefficients? In the course of the 
conjugate gradients procedure (see references), the conjugate directions v* 
of G* (A t D” 1 A + C -1 + Q7 1 ) are generated (Q7 1 is the variance covariance 
matrixof the initial estimates, as explained in paragraph (2.18) part (b)), 
together with the scalars v T * Gvn = a k . The conjugate directions have the 
property 


v i G Vp = 0 


for k ^ p . The estimator Implied by R Iterations of this procedure is 
c * Gf(A T D _1 m + Q" 1 c.) = G A T D" 1 (z + n) + G Q7 1 (c + Ac.) 


where c, is the vector of Initial estimates, Ac, the errors in this vector, and 

* . 
r* f _ I T 

G = 2_ Vk Vk 
k si 

The variance-covariance matrix of the updated errors, for ordinary least 
squares, is (assuming that E {n Ac, | = 0 ) 

E j(G A T D _1 n) (G A T D" 1 n) T }+ E-f(G Q^Ac.) (GQl 1 Ac, / } - GA T D -1 D D -1 A G 

+ GQl l Q. Q7 1 <3 = G (A T D' 1 A + Q', 1 ) 3 = GGG 

R T R t 

= ra; v^ or k v k a'k 1 » Y a 7 VkV k =G 

k = l k=l ~ 


Therefore, the variances of the errors In the 
ding diagonal elements of G : 



K 



k = l 


3 

v*. b« a 



are equal to the correspon- 


where Vk, Bt a is the element of v corresponding to d Bi in d (v and d 
have the same dimension). The same result applies to least squares col¬ 
location. All that is needed to obtain the , according to the formula 
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above, is knowledge of the a. k and of the v * , which can be saved as they are 
created, during the klh iteration of the procedure. For the efficient algorithm 
to be applicable, Q” 1 must have a suitable structure. One case is when 
Ql 1 is diagonal, or can be satisfactorily approximated by a diagonal matrix. 


2.17 The Error Matrix in the Band-Limited Case 


From (2.41) and (2.77) results 

E t = C - Co, (C„ + Df 1 CL = C - Fcl* - C - (A T I)" 1 A + C 1 f ! A T D" 1 A C 

= [1 - (A T D _1 A + C J f 1 A T D _1 A]C - (A t D' 1 a +C -1 rt(A T D' x A + C-Va T D" 3 A]C 

= (A t D" 1 a + C -1 )' 1 (2.101) 

In the case of the best unbiased estimator C' 1 is not present in the normal 
matrix, so (2.101) becomes 

E t = (A t D _1 A) -1 = G" 1 (2.102) 


which is the well known expression of the error matrix for ordinary least 
squares. 

Because of the block-diagonal structure of the variance-covariance matrix,^ 1 ' 
the estimates of c“ of different orders are uncorrelated. The diagonal 
elements of the variance-covariance matrix are the variances of the estimated 
coefficients total eri'ors (i.e., sampling plus propagated noise), in (2.101). 

Obtaining the variances of the C nB in band-limited collocation is formally 
identical to getting the variances of the estimates in ordinary least squares, 
according to (2.101) and to (2.102) 


2.18 The Use of a priori Information on the Coefficients 

Assume that all coefficients up to some degree and order M are 
approximately known, and that Q, is the variance-covariance matrix of their 
errors. This could be the case where a model of the gravity field has been 
obtained from data gathered using artificial satellites, complete to degree M , 
and terrestrial data is to be used to Improve the existing coefficients and obtain 
new ones beyond degree M . 

Three possible approaches to this question will be discussed here, using 
a "point values" formulation in the first two cases for simplicity. 


When the data set is not sparse. 
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The terrestrial data can be used separately to obtain a model to degree 
and order M T > M , together with the variance-covariance matrix Qr of 
the coefficients' errors. To combine the satellite and the terrestrial coef¬ 
ficients one can set up the following observation equations 






I.* 

! 

+ 


a 


N 1 

ol 


■1. 


.It. 


(2.103) 


where c is the vector potential coefficients and I, is the (M + l) 2 x (M + l) a 
unit matrix augmented with zeroes on the right, and I T is the (Mt + l) 2 x 
(Mr + l) 3 unit matrix, while c, is the vector of batellite'teoefficients and 
c t the vector of'terrestrial' boefflc ients; £ and _t are the corresponding 
vectors of residuals. The best linear unbiased estimator for the combined 
system of observations is 


£ - (til It 1 


fQi 1 o TI.1 f 1 [I, I r ] TQf O Tc.-j 

Lo Q^JLitJ Lo Q7 l IcJ 


- + QtV 1 (Q'.'c. + Q ? 1 C t ) 


(2.104) 


If Q. and Q T are the Inverses of ordinary least squares normal ma¬ 
trices G, and Gt , then the error matrix of the combined solution 

E = (Q".‘+ OTt 1 ) -1 = (G. h-Gt )" 1 

corresponds to the propagated noise only. If they are "collocation" 
matrices of the type (A T D" 1 A + C’^see expression (2. 77)), then the error 
matrix includes the effect of the sampling error as well. Most satellite models are, 
to date , 'least squares-type" and It would be incongrous to combine them 
with "collocation-type" models, terrestrial or otherwise. The problem need 
not be a serious one, because geodetic spacecrafts so far have orbited at 
altitudes of 800 km or more, where the field is much smoother than at the 
surface, so the sampling errors are bound to be small compared to the pro¬ 
pagated data errors reflected by Q, . 

In the case where the terrestrial model has been derived from a reg¬ 
ularly sampled data set using the "band-limited approach" of previous para¬ 
graphs, Qt is block-diagonal, and those blocks corresponding to orders 
m > M ar$ identical to the corresponding blocks in (Q7 1 + ) -1 . From 

this it is not difficult to conclude that the coefficients in the combined model 
up to order M will be somewhat different (and presumably better) than those 
in either the satellite or the terrestrial sets, while those above M will be 
identical to the corresponding terrestrial coefficients. 






(b) A priori values included as data in the adjustment : 
Consider the system of observation equations 


r-i 


r 

>1 

i ■ 

•A-i 


+ 


n 


p 

0*1 


1 

ml 

_i 


.1.. 


(2.105) 


which is the system (2.66) augmented with equations of the type 

^m(S) - e M(a) = C a ? 

where e^(s) is *h e error in the "satellite model" coefficient C. 

The normal matrix of collocation is 

E7 1 * ([A T il] TD" 1 O “IfA-j + C' 1 ) (2.106) 

-o Qw 1 JLi.J 

and the optimal estimator of the band-limited type is 

c = (A r D" 1 A + Q7 1 + C"V 1 (A T D" 1 m + Q? c.) (2.107) 

Once more, if the data in m has been sampled regularly on the sphere, 
the estimated coefficients will be affected by the existence of a priori values 
only if their order is no higher than M . Naturally, this is a desirable sit¬ 
uation. Also It is Important that the error measures corresponding to E T 
and Q« be congruous, though this is probably not very Important in the case 
of satellite models obtained with high-orbiting spacecraft. Notice that (a) 
and (b) are equivalent when Qt = (A t ET 1 A) -1 and when C -1 is excluded 
from (2,106)-(2.107). In other words: these first two approaches are equiv¬ 
alent for ordinary least squares. 

(c) The method of Kaula and Rapp ? 

W. Kaula (1966) proposed a technique for simultaneously filtering er¬ 
rors out of a terrestrial data set and improving the coefficients of a satellite 
model. This method was later developed by R. Rapp (1968), who more recently 
(1978) used it to Improve a global data set of mean 1° x 1° anomalies by com¬ 
bining it with the potential coefficients of the GEM-9 model. This adjusted 
data set was used by the author of this report to create the 5° x 5° mean anom¬ 
alies analysed in one of the numerical experiments of section 3. 

The idea is to satisfy condition equations of the type 

' (4Try (n - 1) j3 a S 9 *") -1 Z L I (8. *)<*<* (Agy + v t j) * 0 

i*o j*o*fr t} (2 . 108 ) 
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while minimizing the quadratic form 


g(d, v) = d T Q - , 1 d + v T D" 1 v { 2.109) 

where v is the vector of corrections v t i to the mean gravity anomalies 
Ag i j , and d is the vector of corrections to the satellite potential coef¬ 
ficients 6 £* b (s); & is the nth degree Pellinen factor discussed in paragraph 
(4.3), y is the mean value of equatorial gravity, and S is the ratio between 
the radius of the Earth's largest inner geocentric sphere and the mean Earth 
radius. In matrix form, the condition equations (2.108) are 

c, + d - A t Ag - A t v = 0 (2.110) 

where A is a 2N 2 x (M + 1) 3 matrix having columns a^„ of the form 

a p “ = (4n(n - 1) yS*® A)" 1 if Y n “(0,A)da. . . f Y B “(9,A)dcr] T 

•)Oq o vO N_1 3M_1 

so, except for the factor (4rr(n - l)yS n '* ;3 jS n ) -1 , A is the "area means version" 
of the matrix of system (2.66), and has the same properties as the "A" mat¬ 
rices considered so far. 


The optimal estimates are given by the expressions 
c = c. + d (2.111-a) 


= Ag + v 


(2, lll-b) 


where 

d = -((A D A)" 1 + Q?)~ l (A t D A)“ l (c, - A T £g) 


(?. 112 ) 


and 

v = I) A Q, d (2.113) 

If the data set is both complete and of uniform quality, the matrix A T DA 
has the block structure first discussed in paragraph (2.15). The presence of 
D instead of D -1 makes no difference to the calculations needed to set up 
and invert the matrix: the procedures are those already explained. Ter¬ 
restrial data, however, is usually both scattered and of varying quality 
(i.e., different noise variances). For this type of data, therefore, the meth¬ 
ods for scattered measurements given in paragraph (2,16) could be used. 


2.19 Optimal Estimation over a Band of Spatial Frequencies 
Assume that the signal is of the type 
m = A c + n 
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where A and c may now he infinite (l.e., all degrees from 0 to * 

maybe present). If c' is a sub-vector of c comprising, say, the first 

1? coefficients, and if £ is a vector of estimates of a function s at 
a given set of points on the sphere, such that the values of s depend only on those 
of c' according to the relationship 

g « Be' (2.114) 

where 8 is some matrix of appropriate size, not necessarily of the same type 
as A , then the optimal estimator for £ is, according to (2.39), (2.75), 

i - C,.(C,* + Df 1 m 

= M{s z T l (On + Df 1 m = M{B o' c' A (C„ + Df x } 

= Be'A'’(C M + D)" 1 m (wnere C’ = M-fc' c ,T l y 

= B C 0 - x (C ZM + Df 1 m L J 

or 

| » Be- (2.115) 

according to (2.40). 

Expression (2.115) indicates that the optimal estimates of a band- 
limited function s from data m are identical to the values of a obtained 
from the optimal estimates of the coefficients c' by means of the rela¬ 
tionship 1= b£ . 
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3. Numerical Examples 

This section presents several computed examples to Illustrate some 
of the ideas and methods discussed earlier on. The question posed here is, 
basically, that of the accuracy of the various procedures, and Is answered 
by means of error analysis carried out with the formulas for error variance 
developed in section 2, and also by analysing simulated data and comparing 
the recovered coefficients to the original ones in order to find the actual 
errors. A comparison of the rms of these errors with the theoretical rms 
(i.e., the square root of the variance) provides both a check on each set of 
results and, more important, shows just how adequate an error measure 
the theoretical rms can be. Besides error analysis and simulations,this 
section shows, in the last paragraph, the results of the harmonic analysis 
of a real data set: a 5° x 5° equal angular set of mean gravity anomalies 
covering the whole Earth, from which the coefficients of the disturbing 
potential have been recovered to degree and order 36 by means of least squares 
collocation, using the "Toeplitz matrix" approach of paragraphs (2.10) and 
( 2 . 11 ). 


3.1 Generation and Analysis of Simulated Data 

As explained in the proceeding section, the variance of the error in 
the estimate of C® depends on the power spectrum (or covariance function) 
of the signal, and on the variance-covariance matrix of the noise. The pro¬ 
pagation of the noise is quite straightforward, and anybody who has had any 
practical experience with adjustments of geodetic networks and the like al¬ 
ready has enough "feeling" for this part of the error measure, and is cap¬ 
able of understanding its significance when its value is given to him. The 
part corresponding to the sampling error is somewhat different; it involves 
a rather unusual geometric average over rotations, and this type of error 
measure, while not exactly new (collocation, based on this measure, has 
been around since the mid-sixties) is not so familiar to geodesists yet, and 
its use,in harmonic analysis In particular, far from common practice. For 
this reason, it is probably fair to the readers to provide some illustration of 
how "close" this part of the error measure is to the actual sampling error 
that occurs when data of the assumed power spectrum is analyzed in any of 
the ways discussed so far to recover spherical harmonic coefficients. By 
"close" one means that the actual numbers measuring the theoretical and the 
actual variances (or rms) should differ from each other by a small percentage, 
or some equally clear-cut criterion. 

The theoretical variance considered here is the variance of the esti¬ 
mation errors per degree o^ B , defined in terms of the error measure of 
section 2 as follows (see paragraph (2.8)) 






n 

- I 

u = 0 


I a ) 

l I 


a—o ■ = o ot-° 


(2n + 1) 


^£.a*tt >z _f a8 + (fpaJtCjj + H)_fBa) 

(3.1) 


The rms of the error is the square root of this variance, and the ratio of 
this rms to the rms per coefficient: 


- <*- 




l l 

«=o a=o 


(2ia,«i®- &?>’ (C 


SZ 


+ D) f, 


a p 
^m))® 


(3.2) 


multiplied by 100 , or the percentage rms error per degree, is the theo¬ 
retical quantity to be compared to the "actual” percentage rms error per 
degree derived from the analysis of simulated data with the same statistical 
characteristics (i.e., a a, » , C, t , and D) in formula (3.2). 

The cr a „ were computed with subroutine NORMAX (Appendix B). 


To obtain the actual percentage rms error per degree, sets of simulated 
data were created on full regular grids as follows! the artificial data con¬ 
sisted of area means computed globally, using the algorithm outlined in para¬ 
graph (1.7) and subroutine SSYNTH (Appendix B), on the basis of expression 
(1.2). The C“ , complete to degree and order N„ x ^ N = came from 

sequences of random numbers. The random numbers, obtained using the IMSL 
subroutine "GGNOR" with generatmg'seeds"of the order of 10 4 , were scaled 
to give them the desired degree variances cr B a . For each simulation, a se¬ 
quence of (N BaX + l) a numbers was obtained, the first corresponding to c 0 o » 
the second, third, and fourth to C 1 0 , C u , and S n , respectively, and so 
forth. If r n0 , r® x , . . . r n ® were the (2n + 1) numbers corresponding to 
degree n , then the scaling that resulted In the corresponding C n a , was 



(3.3) 


The harmonic coefficients C® obtained in this way were the ''actual" coef¬ 
ficients to which the fc®, recovered by some of the procedures described 
in section 2,were then compared to obtain the actual percentage rms errors 
per degree 

C» = l t I (££ " x 100 (3.4) 


The analysis of the simulated data was done with subroutine HARM IN (appendix 
B). This type of numerical experiment was carried out three or more times 
in each case, varying only the seed used to generate the random numbers,much 
as a Montecarlo-type of analysis is conducted. The seeds were chosen widely 
apart, to ensure that the correlation between "trials” would be virtually nil. 
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The maximum degree and order in the set of artificial coefficients, 

N Ba x , was chosen so that the power in the mean values above degree N Ba » 

8000 

Pn..x = £ ft V? (3.5) 

k>tl 

were less than 1% of the power between degrees 0 and 2000. The 
are the Pellinen coefficients»discussed in paragraph (4.3), corresponding 
to 5° x 5° area means. The values used for the degree variances <T n had been 
empirically obtained from terrestrial data in the manner described below. 

The data were supposed to be noise-free, as only the sampling part of the 
error was studied in this way, for the x'easons given at the beginning of this 
paragraph. The estimators being linear, the propagated noise and the sampling 
error merely add arithmetically to each other, and can therefore be studied 
separately, if so desired. Summing up, it can be said that the simulated 
data consisted in global data sets of artificial gravity anomalies, averaged 
over equal angular grids. 


The empirical degree variances were obtained as follows: up to degree 
100 they were those Implied by a set of coefficients, complete to degree 180, 
obtained by R. Rapp and associates at O.S.U. from a global data set of 
1° x 1° mean anomalies. Above degree 100, the (Ag) were obtained from 
a model of the form 


(3.6) 


(Moritz, 1976), where the parameters a 2 , a s , S a , S a , A and B have 
been adjusted to fit existing gravimetric data, satellite altimetry, satellite 
field models, and other geophysical data. The parameters used in most 
examples were 


0c x = 3.4050 = 0.998006 A = 1. 

a 3 = 140.03 S 9 = 0.914232 B = 2. 


corresponding to the best model of this type given in a report by R. Rapp 
(1979) who, in the same work, discusses also the empirical degree variances 
obtained from his 180, 180 field model. The degree variances implied by (3-6) 
with the parameters listed above are also very similar to those obtained by quite 
different means by Wagner and Colombo (1979), who analyzed the (Fourier) power 
spectrum of short arcs of GEOS-3 altimetry, and converted their average to 
a spherical harmonics spectrum using formulas that follow from the 
relationship between spherical harmonics and Fourier series. The empir¬ 
ical variances for n £ 100 are included in the listing of subroutine NVAR, 
in Appendix B. 

In order to understand how critical the choise of empirical degree variances 
is to the theoretical and actual errors, a different two-term model obtained 
by C. Jekell (1978) was used as well. This model has the following parameter 


- 70 - 






values 


= 18.3906 Sj = 0.9943667 A = 140. 

Ola = 658.6132 S 3 = 0.9048949 B = 10. 

All the examples considered In this section refer to complete equal angular 
sets of mean values with the same statistical properties of terrestrial gravity 
anomalies (as far as such properties are known). The analysis of actual 
mean gravity anomalies is shown in the last paragraph. Besides being im¬ 
portant in geodetic studies, gravity anomalies constitute a type of geophysical 
data with reasonably well known statistical properties, and their study here 
is meant to give the reader some idea of how effective are the ideas presented 
earlier when it comes to handling "real data'' ( or something resembling it). 


3.2 Agreement between the Actual and the Theoretical Measures of the 
Sampling Errors 

Table (3.1) lists side by side the theoretical precentage rms per degree 
of the sampling error according to(3.2) and the actual value of this percen¬ 
tage for two different sets of coefficients (i.e., from random sequences with 
different seeds). The coefficients were recovered using the quadratures for¬ 
mula 

c?. - 4 -A- E E *»] 

P* 15o 1=0 •'Olj 

This type of formulas has been discussed in paragraph (2.6). The simulated 
data consisted in full sets of 5° x 5° mean anomalies obtained from harmonic 
coefficients complete to degree and order N»»* = 140 . The power above 
degree 140 in 5°x 5° anomalies is negligible, according to the empirical 
power spectrum model that was used. The results shown here are fairly 
typical of similar tests conducted with other quadrature formulas, so the 
conclusions that can be drawn are likely to be valid for the analysis of area 
t means by numerical quadratures in general. There is clear agreement 

| between the theoretical and the actual rms of the errors, and not just the 

[ average rms of actual errors, but the actual rms of each trial as well. The 

agreement is close, and the reader will probably agree that to use a theoretical 
error measure that can predict the actual error so well is a meaningful way 
of quantifying the error. 


3.3 Accuracies of Various Quadratures Formulas 

Five quadratures formulas for area means have been studied: the first 
four of the type 

x~i aw-1 


h-i aw-L * _ 

Z I Y® (9,A)dq Ag„ 

i = o :=oJ<7 u 
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Table 3,1 

Comparison of actual versus theoretical percentage rms 
error per degree. 

5° x 5° mean anomalies, N„ x = 140, 0. mgal rms noise. 


n Actual, No 1 Actual, No 2 Average Theoretical (expres- 
(seed=53218) (seed = 31765) 1 and 2 sion (3.2) x 100) 


2 

! 0.50 

0. 62 

I 0.56 

0.51 

5 

1.06 

1.10 

! 1.08 

1.23 

10 

4.99 

5. 81 

5.40 

4.89 

15 

13.06 

11.03 

12.05 

12.75 

20 

23.48 

25.97 

24. 72 

24. 77 

25 

23.04 

30.82 

26.93 

30.03 

30 

45.84 

46.26 

46.05 

43.28 

*36(N) 

55.33 

61.10 

58.22 

60.35 

40 

70.61 

75.82 

73.22 

75.14 

45 

85.91 

86.41 

86.16 

87.03 

50 

101.51 

101.43 

101.47 

101.47 


differ only in n n : 

(a) /x n = the optimal de-smoothing factor given by expression (2.3G) in para¬ 

graph (2.7): 

(b) Mn = 

(C) xrfi; 

(d) 

(e) The optimal quadratures - type formula, in the least 

squares collocation sense (i.e., minimum combined error measure) for 
the given grid; si<” >1 and noise (paragraph (2.8)). The grid was equal 
angular in all five ,es. Table (3.2) compares the percentage rms of the 
errors per degree for a 30° x 30° grid; table (3.3) corresponds to a 10° x 10° 
grid ; and table (3.4) to a 5° x 5° grid. N„, x was 100 for the first two tables, 

'*n«! 140 for the last. All these values are theoretical, computed in accordance 
r *V formulas in paragraph (3.1). In all three cases noise is not present, 

• ■ • rr >rs are purely sampling errors. 

• i : i 'i. :>, corresponds to a 5° x 5° grid and an uniform noise of 
'■ *!-«• • ffect of the noise has been included in the results. 

. < yoe correlation coefficients for the various methods, 

• . ■ ■ "'relation coefficient for the nth degx’ce is defined as 


. »• 


a «d. 











(3.7) 


■ = o = o 


p.-it t cf. 4“ (t r<£v<£ 

« = ° QpO * = ° OF 

and It is also equal to 

f 

Jo Ag„ Ak.dcr 


A = 


/a ^ dff /a 


(3.8) 


This coefficient can be regarded either a a measure of the agreement between 
the actual and the recovered coefficients of the nth harmonic, or of between 
the nth harmonic Ag„ in the signal and the harmonic Ag B that can be computed 
from the recovered coefficients, if the C* could be seen as random variables 
with gaussian distribution, the interpretation of ft would move along well-worn 
paths; however, as it was mentioned in paragraph (2.5), there are some unex¬ 
pected problems when extending the idea of a gaussian random process to the 
sphere, so is better to chose another approach. One could regard the coef¬ 
ficients of the nth harmonic as the coordinates of a vector in (2n + 1)-dimen¬ 
sional euclidean space. The actual coefficients will define thus one vector, 
and the recovered coefficients another. Expression (3.7) then merely defines 
ft, as the scalar product of this two vectors. Likewise, expression (3.8) is 
that of the scalar product of two elements of a function space. The angle 
formed by these vectors is 0° when correlation is (maximum), and 90° for 
0 correlation; a minimum correlation of -1 corresponds to the case when 
the vectors are equal but of opposite sense. The scalar product is independent 
of any scale factors that may multiply the vectors: it depends only on their 
mutual orientation. For this reason, the correlation coefficients are the same 
for the four quadratures formulas of the type 


« 


N*1 

1 = 0 


l 


SN- 1 

r Agn / C (e,X)dc 

J = o 


because the difference between the errors for the same harmonic, predicted with 
two different formulas of this kind, consists in a scale factor MavMc • 

Clearly, as the rms of the error Increases with n , ft, decreases from almost 
1 where the error is smallest (very low degrees),to below 0.5 where the 
error exceeds 90% (highest degrees analyzed). 


Observing the percentual rms of the errors in the first four tables, it 
is easy to see that they by no means reach 100% as soon as the Nyquist 
frequency is reached (n = N), but that they remain substantially smaller than 
100% even at degrees considerably higher than N ; this is in line with the 
conclusions in paragraph (1.3). The optimal estimator itself cannot have an 
error larger than 100% , be it due to sampling, noise, or both. Otherwise, 
a null estimator (one that predicts only zeroes) would be better than the op¬ 
timal, which is not possible. 


Table (3.7), compares the theoretical errors with zero noise (i.e., 
the sampling errors) of the collocation estimator obtained, first according to 


- 73 - 






< 7 a 3 implied by R. Rapp's model (used in all the other tables), and then ac¬ 
cording to the "2L" model of C. Jeleky, both described in paragraph (3.1). 

This table is included here to give the reader an idea of how sensitive the theo¬ 
retical error variances are to the empirical degree variances used to compute 
them. The '’2L" model has considerably more power than Rapp's at high degrees, 
and this may be reflected in the somewhat larger errors in the corresponding 
col umn of the table. 


Table 3.2 

Theoretical percentage rms error per degree. 


30° x 30° mean anomalies. 


= 100, 0. mgal rms noise. 


Optimal 

Estimator 



13. 

,17 

25. 

,60 

65. 

,13 

96. 

,01 

99. 

,27 


_1_ _L_ _I_ 

4 TT 3n 4n infi 


13.57 
28.72 
101.41 

194.08 97.13 882.98 

915.49 105.44 17835.58 



Table 3. 3 

Theoretical percentage rms error per d egree. 

10° x 10° mean anomalies, N ltx = 100 , 0. mgal rms noise. 


Optimal 

Estim. 


1.52 
2.24 
4.27 
8.74 
13.19 
31.60 
40.43 
41.41 




* 

Mb 

4tt Pn 

1 

4 TT 

1 

4rr~8f 

1.63 

1.71 

2.04 

1.65 

2.46 

3.00 

4.76 

2.56 

4.68 

5.83 

9.51 

4.94 

9.60 

10.99 

16.39 

10.16 

14.02 

15. 89 

23.68 

15.19 

33.03 

33.10 

37.84 

37.43 

41.41 

41.41 

46.45 

49.91 

42.18 

42.36 

51.41 

53.12 

59.94 

61.26 

63.41 

88.17 

78.99 

94.69 

79.07 

158.26 

78.54 

89.75 

79.05 

167. 49 

87.40 

116.57 

87.59 

258.09 

94.03 

164.50 

97.85 

443.76 

95.98 

193.17 

100.35 

646.88 

97.30 

232.24 

101.88 

1010.13 

98.30 

303. 89 

103.49 

1828.98 

97.40 

281.73 

97.77 

2619.98 

98.08 

485.72 

98.45 

' 8694.11 















Table 3.4 

Theoretical percentage rms error per degree 


5° x 5° mean anomalies, 


Optimal 

n Estim. u. 


140 , 0. mgal rms noise. 


0.47 

1.36 

9.74 

15.14 

30.58 

42.88 


0.49 
1.44 
10.11 
15. 55 
31.06 
43.28 
59.15 
73.51 
83.62 


Table 3.5 

Theoretical percentage rms per degree 

5° x 5° mean anomalies, N a » x = 140 , 5. mgal rms noise. 


n 

Optimal 

Estim. 

2 

8.74 

6 

9.33 

12 

32.80 

18 

35.94 

24 

53.32 

30 

61.97 


Table 3. 6 

Correlation factor per degree. 

5° x 5° mean anomalies, N ns ,= 140, 0. mgal noise. 


n Optimal Estim. Simple quadratures 


2 

6 

12 

18 

24 

30 

36 (N) 
42 

48 
















Table 3. 7 

Comparison between the theoretical errors (optimal esti- 
mator)when different covariance functions are assumed. 
N„* =* 140 , 5° x 5° mean anomalies, 0. mgal noise. 



The results in the tables show that the optimal estimator errors are the 
smallest, as expected, and that the quadratures formula with the optimal de¬ 
smoothing factors is the best of the four simple quadratures formulas compared here, 
though it is not quite as good as the optimal estimator, also as expected. In 
the case of zero noise each of the three non-optlmal quadratures formulas 
has errors that, in some region of the spectrum, are smaller than those for 
the other two; concretely, the de-smoothing factor works better 

for low degree harmonics (or n s; -JrN> i.e., the percentage tins of the 
error is less than for the other two formulas there, the factor p a = 
is best for middle harmonics (approximately l/3N<n £ N), and 
/i B = is best above the Nyquist frequency N . It is possible, therefore, 
to obtain a simple "composite" quadratures formula that combines the good 
properties of all three formulas, by defining its de-smoothing factor as follows; 

| )3 a if 0 * n * 1/3N 

Ms = r^T w ^ere 7j a * i /3 a if 1/3N <n s N (3.9) 

Vb (l if n > N 


This composite formula has been implemented in the version of-subroutine 
HARMIN listed in Appendix B, through the subroutine can be easily changed 
to compute other quadrature formulas. 


It is clear from the tables that, while better than all the others, the 
optimal formula is only marginally so; from a practical point of view, the 
simple "composite" formula (3.9) above Is virtually as good, but it is much 
easier to implement and compute. Therefore, when analysing data of the type 
considered here (resembling mean anomalies with uniform noise on an equal 
angular grid), the quadratures formulas discussed the composite In 
particular, are about as good as any linear technique for estimating the coef¬ 
ficients, and also very easy to program and very efficient. 








There may be cases, however, when this is not true. If the data set 
were very noisy and/or Incomplete, or if the coefficients to be recovered were 
not those of the signal, but those of some complicated transformation of it 
(as in the case of satellite-to-satellite tracking data, where the signal depends 
on a combination of differences ox raaiai and horizontal derivatives of the grav¬ 
itational potential) so no simple Integral formula like (1.4) exists that can be 
readily discretized into summations like (1.5) or (1.6), then simple quadrature 
formulas like those in section 1 would be of no real use, and the optimal es¬ 
timator could provide the only practical way of obtaining (he coefficients. Both 
collocation and least squares adjustment are very similar, as shown in section 
2, and both can be implemented with reasonable efficiency if nothing simpler 
is available. More important, the optimal formulas provide a theoretical 
background on which one can build a coherent and comprehensive understanding 
of the other linear techniques for spherical harmonic analysis. It is, after all, 
because of the theory developed in the previcu s section that it has been possible 
to obtain the results shown in the proceeding tables, results that constitute the 
factual basis for these considerations. 

The techniques for setting up and inserting the variance-covariance 
matrix of the data are of interest In a number of estimation and filtering 
problems, besides harmonic analysis. The author hopes that the examples 
provided in this section will encourage the wider use of least squares col¬ 
location for the processing of large, global sets of data, both point values 
and area means. Creating the simulated mean anomalies with subroutine 
SSYNTH, obtaining the "weights" x? of the optimal quadratures formula with 
NORMAL^ (Appendix B), and recovering the b®, for comparison with the or¬ 
iginal C„ up to degree and order 72 (same 5000 coefficients) on a 5° x 5° degree 
grid (about 2600 "data" values) took less than 20 seconds, central processor 
unit time, in the AMDHAL computer at Ohio State. To recover the harmonic 
coefficients up to degree and order 130 from a complete equal angular set of 
1° x 1° mean values would require less than two hours, using the same machine. 
Most of this time would be dedicated to creating and inverting (C tl + D). 

Butin the optimal estimation and filtering of geoidal undulations, deflections 
of the vertical, and any other function of the gravity field estimable from (say) 
gravity anomalies, the fact that (C,* + D) can be set up and inverted efficiently 
transcends harmonic analysis. A major Implication Is that such extremely 
large global adjustments require a computational effort that is already within 
the reach of most researchers. 

As mentioned, already In section 2, matrix (C lf + D) may become poorer 
conditioned, i. e., its numerical inversion less stable, as the data distribution 
becomes denser. This tendency towards instability was noticed: when there was 
no noise (D = 0) the R(m) matrices had to be regularized by adding a small 
positive constant k to each diagonal term (Colombo, 1979, par. (4.5)) before 
they could be inverted successfully. This constant, which In most cases was 
much smaller than the diagonal elements It was being added to, was lCf 6 for 
30° x 30° and 10° x 10° , but had to be increased to 10"° in the case of a 
5 C x 5° grid. When noise was present, the nonzero diagonal elements in D 
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were sufficient to provide stability, and no regularization was needed. Be¬ 
cause of this tendency to instability, the "band-limited" approach of paragraphs 
(2. 13) to (2. 17) may be preferable, whenever it can be properly applied. 

As shown in Table (3.1), the theoretical and the actual sampling errors 
are almost the same in most cases. The propagated noise measure is very 
easy to compute in the case of imcorrelated noise, and can be added to the 
actual sampling error (variance) to get an estimate of the total error, both 
actual and theoretical. This estimate, where the sampling part is the result of 
a Montecarlo-like approach, is much easier to obtain than the theoretical one 
that involves setting up matrix (C 2Z + D), or at least the R(m) matrices. In the 
case of equal angular data sets like those considered here, this empirical es¬ 
timate is likely to be just as accurate, when it comes to judge the performance 
of any given type of harmonic analysis. Such estimate has been used to evaluate 
the likely errors in the potential coefficients obtainedJ^rom 1° x 1° mean anom¬ 
alies using the quadratures formula £ 0 /o u (9»X)da /Sgij . 

The Montecarlo method described in paragraph r (5". 1) was implemented with the 
help of subroutines SSYNTH and HARMIN, the error variance being equated 
to its average over three "trials" (three sets of coefficients created from dif¬ 
ferent random sequences), C. Jeleli, also at O.S. U., who undertook this work 
as part of his own research.fitted a quartic to the percentage rms per degree 
thus obtained. When this was expressed as a function of the "normalized de¬ 
gree” n/N , the quartic fitted equally well the theoretical results for 30° x 30° , 
10° x 10° , and 5° x 5° presented in this section, jekeli's quartic expression 
for the truncation error is (private communication): 

-^5 x 100 = [(-16.19570 (-J) + 30.34506) (-^) + 40.29588 ) (^-) 2 
n (3.10) 

It is quite remarkable that such a complex phenomenon can be described satis¬ 
factorily by such a simple law. 

The expansion of the simulated 1° x 1° mean anomalies was complete 
up to degree and order N,» x = 300 . Creating (or analysing) area mean values 
up to degree and order N BBX = 300 required about 50 seconds o.p.u. time 
using double precision arithmetic in the AMDHAL computer at O.S. U. 

Table 3.8 shows the actual percentage rms sampling error per degree 
as computed in one of the trials, and the percentage rms propagated noise 
(theoretical) corresponding to a 1 jngal rms noise in the data. Clearly, the 
errors are much smaller than for any of the cases considered previously; 
this improvement is due to the finer sampling (the sampling error tends to 
zero as the area of the blocks tends to zero). The data, however, tends to 
be noisier when averaged on smaller blocks, so the propagafed noise may in¬ 
crease. For a given rms error In the data, multiply the number in the "pro¬ 
pagated noise" column by this rms (in ragals) to obtain the corresponding per¬ 
centage. These numbers arc only valid for the estimator where p n = —1— . 
Repeated trials with different random coefficients resulted in much 4tt £» 
the same percentages for the sampling errors, so these values are probably fairly 
typical. -78- 
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Table 3.8 

Percentage rms errors per degree for a 1° x 1° 
equal angular grid. n n ■ , N a , x = 300, 

€ Ag = 1 mgal 


actual 

propagated 

n 

sampling 

error (theor.) 

2 

0.01 

0.47 

10 

0.10 

0.72 

20 

0.46 

1.39 

30 

0.98 

1.72 

40 

1.55 

1.93 

50 

2.68 

2.07 

60 

4.71 

2.20 

70 

6.94 

2.34 

80 

9.43 

2.47 

90 

13.51 

2.62 

100 

14.79 

2.77 

110 

17.98 

2.96 

120 

25.06 

3.08 

130 

25.33 

3.26 

140 

30.87 

3.43 

150 

37.02 

3.62 

160 

44.17 

3.81 

170 

44.43 

4.03 

180 (N) 

53.07 

4.24 

190 

61.87 

4.47 

200 

65.90 

4.73 


A possible way of bringing the sampling error down is to use weighted 
area means of the form 



where Ag{ k ) is the kth measurement inside the block a tJ , and where the 
W k are functions of the distance of Ag® to the center of the block. It the 
W k decay gently towards the border of the area element, the resulting 
weighted means will be smoother than the ordinary area means considered so 
far (all W k = N^j) and their harmonic content above the Nyquist frequency 
will be atenuated. Consequently, the harmonics below N can be recovered 
with less sampling error. This Idea certainly deserves further study, ob- 

vkrnsly, it is applicable only in those oases where the original measurements 
Aare available. 


L 
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3.4 The Analysts of a Global Data Set of 5° x 5° Mean Anomalies 


As a final demonstration of the use of optimal estimators, this paragraph 
presents some results obtained by analysing a real data set, consisting of 5° x 5° 
mean gravity anomalies. These anomalies were obtained from a global set 
of 1° x 1° mean values created by R. Rapp and associates (Rapp, 1978) from 
the combination of land and gravity measurements, satellite altimetry over 
the oceans, and the GEM-9 satellite model (Lerch et al., 1979). 

The 5° x 5° values, £g(s°) » were obtained from the 1° x l e values, 
figtf), using the formula 

= *25" £ **(*)*« 

where I represents summation over all 1° x 1° blocks inside a 5° x 5°block. 
The variances of the £g^ were obtained from the formula 

t,a Vf <skf £ 

where the assumption has been made that the errors in the 1° x 1° values 
were uncorrelated, which is not quite true, as the values used ware the pro¬ 
duct of an adjustment. The variances 0 3 Ag.gO. were different from 5° block 

to 5° block, but they were " homogeneized" as described in paragraph 
(2.9), In order to obtain a (C„, + D) matrix that was easily invertible and re¬ 
sulted in a quadratures-type estimator inexpensive to implement. 

Figure (3.1) shows a comparison between the collocation solution, com¬ 
plete to degree and order 36 , and the same coefficients obtained by R. Rapp 
by numerical quadratur es from t he original 1° x 1° values. The figure shows 
the rms per degree ( y<r,?/(2n + l)) for potential coefficients (cr n a (T) = 
of (Ag) y 3 (n-1)"® , y = 979800mgal). The circles correspond to Rapp's re¬ 
sults, and the triangles to those obtained using collocation as already explained. 
Because the grid used by Rapp is much finer, the corresponding results are 
likely to be less affected by sampling errors, at least up to degree 36 , than 
those obtained from the 5° x 5° anomalies; for this reason Rapp's rms values 
are regarded here as the "true" ones. The solid line corresponds to c. Jeleki's 
"2L" model for the of , used here to obtain the optimal estimator, it is inter¬ 
esting to notice that the triangles follow the circles (or "true" values) rather 
than the line. A common concern among those using this type of estimators is 
to what extent the "a priori" power spectrum or covariance function used to set 
up the estimator may "bias" the results by forcing the spectrum of the output 
to resemble the "a priori" spectrum. Here there is little evidence of such a 
"bias". 


In addition to Rapp's coefficients, those of GEM-10B (Lerch, 1980) 
were also used for comparison. The "collocation" values follow them very closely 
too (the three sets of results agree, in fact, very well with each other). The 
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GEM values are not shown In the figure, to have a cleaner picture. In any 
event, the ’'collocation” values compare very well with the other two, and in the 
higher degrees (31 to 36) where the divergence between Rapp's model and col¬ 
location is largest, GEM-10B fits right in between them. The results shown 
here were presented in a previous paper by the author (Colombo, 197&-b). 


Figure 3.1 
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4. Covariances Between Area Means 


In what follows, the expression "point" covariance refers to 

as 

cov(u(P),v(P')) = E cY Pb(<M * M{u(P)»v(P')} (4.1) 

The covariance between the area means of two functions u and v is related 
to the "point" covariance function by the following integral relationship (Sjoberg, 
1978): 

cov(u u ,v kl ) = t \ f da J cov(u(P),v(P'))da' (4.2) 

where P belongs to the area block a l} , and P' to o kl , while 

ui, = f u(P) do (4.3) 

J Ou 

is the average of u over the block o u of area 

A u = AX(cos(0!) - cos(9i+A0 1 )) (4.4) 

A0 t , A A are the colatitude and longitude spans of the blocks in the "row" between 
colatitudes 0! and 0 t +A0 4 . To simplify the discussion, A0 t =A0 is assumed 
constant; extension to the more general case where A0 t varies from row to row 
is straightforward. 


4. l. Derivation of an Approximate Formula for the Covariance 1 

Expression (4. 2) can be computed by numerical quadratures (Rapp, 1977). 
Because of the double integral, this is a very laborious process, and it is not 
practical if many covariance values have to be found, as in the case of large data 
sets. A more efficient alternative is needed. 

Replacing the covariance function in the integrand of (4.2) with its Legendre 
expansion 

cov(n 1Jf v kl) - da Z q c u /p.<*„) (4.5) 

Here ^ = cos -1 (cos 0 cos0' + sin0 8ln0' cos (A - A')) is the spherical dis¬ 
tance between points P *(0,A), and P' ■(•'.A') in the unit sphere, while 

i ■Bi ou ♦ si §:. 

■ ■0 

is the nth degree variance of the cross-spectrum of u and v. As shown in the 
Appendix, the order of summation and Integration can be reversed. 

1 In this section, expression P„ ( 4>»<) is shorthand notation for P n (cos 
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COV(u,j,V kl ) = —— E Cn" f da f Pn(M d<7 ' 

a=o J a u J a)cI 


(4. 6 ) 


Let P 01 (cob &) be the folly normalized associate Legendre function of the first 
kind of degree n and order m. Applying the Summation Theorem for such 
functions, the P a ( &f> f ) can be replaced in the integrand as shown below. 


cov 


00 u f v * p a 

(UiJ.Vk!) — E L dcr I E P*(cos0) P B .(cos 8 ') cos m(A-A')do' 

a = 0 2 n+l 1 = 0 

(4. 7) 


The sum in the general term of the series above has a finite number of terms 
(n + 1 ), so it is valid to exchange summation and integration once more: 


cov | 


r (5~ij,V| tl ) = E £ L da L Pn«(cos 9) Pa*(cos 8 ') cos m(A - X') d o' 

&ljU kl i:0 Zn + 1 “w k i 

(4. 8 ) 

Writing out the area integrals as colatitude and longitude integrals-. 

, 00 v R 

E -srti E J 

+ 1 »=o . 

•A 3 + AA |*Aj + AA 


cov(u 13 ,v kl ) = 


Al 3 A kl n = 0 2n 

e k +Ae k ^ 


0i + A9i_ _ _ . 

P a ,(cos 8 ) sin 8 d0 
Si 


J 9k+ 9 k P B ,(cos 8 ») sin 8 ’ d 8 > J 

8jc "Aj 


dA 


I 


'X 2 


(4. 9) 
cosm(A-A') dA' 


Calling Aj = A 13 , as all blocks in the same "row" have the same area, we 
introduce two functions. 


and 


_ f AA a if m = 0 

m t(2/m a )(l-cosmAA) otherwise, 

w»(^ T ) i ir" v,cM9,stosde 


(4.10-a) 


(4.10-b) 


Then expression (4.9) can be written, after integrating and reordering terms, as 
follows: 


S2, oo 

cov(u u ,v kl ) = L Q F(m) E L .,1 In.,k cosm(ArAj) 


(4. ll-a) 


This regrouping of the series is valid because, as shown in the Appendix, the series 
is absolutely convergent. The integrals in the definition of l„, tI (expression (4.10-b)) 
can be calculated very accurately and efficiently with recursive formulae obtained by 
M. K. Paul (1979). Regarding j as being fixed, the last expression is also that of the 
Fouiicr series of cosines of cov(Ui 3 , v kl ), with amplitudes 


a,’ - F( m) £ Im,i I Bljk 


(4. ll-b) 


and phases < 0 , - -mAj . 


f: 


L 
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While (4.11-a, b) are valid for all values of X 3 in the interval 0 £ Xj s 2 rr, 
this egression has to be calculated only for Xj - 0,AX , 2AX,... ,jAX. According 
to the sampling theorem for ordinary Fourier Series, at such regularly spaced X } 
expression (4.11-a) takes precisely the same values as the finite sum 

N 

cov(u w ,v itl ) = a^ cos m( X; - Xi) (4.12) 

1*0 

where N = n/AX 

a^ 1 = E aaNh+« + E a aNh-» (4.13) 

h=0 h =1 

To calculate the covariances, (4.11-a, b) and (4.13) have to be truncated, excluding 
all harmonics above some degree N m ax; (4.12) becomes: 

N * N a»X 

COV(u 1Jt V w ) - COV(u u ,V kl )„ - ETE E (In^Nh + «,l)(I»^N h+ ., k ) 

w max isoL»ao t =i 

K N **x 

+ E E ( I n ,aNh-m,l)( ln, 3 Nh-»,k)l F( m ) cos m(Xj - Xj) (4.14) 

h= 1 n:i J 

where (2K +1)N s N mflY . 


4.2. Choosing N 


max 


To calculate the values of covariances between equal angular mean anomalies, Ag, 
the value of N may could be chosen so that the percentage error 


v = 


COVfSgtj.SgkOq - COV(5g 1J ,I& 1 )N 


C0V(2g u t £g kl ) q 


2XJL 


x 100 


(cov(Ag M ,Ag k i )q being computed by numerical quadratures) does not exceed a pre¬ 
scribed upper bound e . The smallest N rriny that meets this condition increases 
towards the poles, because the decrease in area of equal angular blocks with latitude 
means that the averages have a high frequency content that increases accordingly. 

On the other hand, the absolute values of the integrals of the P u , and therefore 
the Fourier coefficients a, in (4.12), decrease quite fast with increasing order m 
near the poles, so their contribution soon becomes insignificant. This is fortunate, 
because the need for lengthier calculations for each Fourier coefficient nearer to 
the poles can be offset by the existence of fewer coefficients there. In fact, oecause 
of the finite arithmetic of digital computers, all a, for blocks less than 30° from 
the poles are rounded off to zero for m considerably less than N in the cases 
presented here. Because of this, calculations near the poles can be less laborious 
than close to the equator, in spite of the larger N mav . To take advantage of this 
in the programming of (4.14), the a, coefficients were not computed above the first 
m for which the following condition was met: 
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a 2 


n-1 AS 

.£ H 

A = o 


* 6 


and 


£ k\ 

{,= 0 * 


6, 


where 6 £ lO" 12 . This ensured no change in the first six significant figures of 
the result when compared to the case where no coefficients were ignored, but 
there were great savings in computing. 


4.3. Numerical Examples 

To verify the accuracy of the truncated series in expression (4.14), mean 
gravity anomaly covariances were computed both with this formula and by numeri¬ 
cal quadratures. All calculations were carried out in double precision (32 bytes), 
some of those for expression (4.14) being repeated in extended precision (64 bytes). 
The agreement between both sets of results was better than 6 significant figures, 
suggesting that both expression (4.14) and Paul's recursives (used to obtain the 
I nBjI ) are quite stable numerically, at all latitudes. 

The covariances computed were between a mean anomaly in a fixed block 
and all the anomalies in the same row of blocks (i. e., all blocks bounded by the 
same parallels), calculations being done for several rows, at close intervals from 
equator to pole. Results for only a few of those will be shown here, because they 
are typical of the rest. Both a 5° and a 1° grid (equal angular) were studied. The 
results were compared to those obtained by numerical quadratures, and by the 
approximate formula 

_ N »«* 

cov (Agij,Ag kl ) =- ^Ai&kCn’ P n (tl>YY<) (4.35) 

n = O 

where Y and Y' are the center points of and while 

A>,i = -- — r — [ P n -i (cos 0 O ) - P n +1 (cos 0 O ) ] 

’ 1 - cos yo 2 n +1 

is the nth degree Pellinen smoothing factor (the formula is Meissl's), and 

Ipo = cos' 1 1" ( COS ( 0J + 1 ) - COS 01 ) + 

This formula gives the covariances of averages on circular blocks of the same 
area as that of the equal angular blocks in the row between colatitudes 8 t and 0 )+1 . 
These covariances are used sometimes as approximations to the equal angular co- 
variances between mean values. 

The numerical quadratures technique consists in the following; (a) sub- 
do .uiug each block with a grid of k equally spaced latitude and k equally spaced 
longitude lines; (b) computing the covariances between point gravity anomalies 
at the nodes of each subdivision (there are k 4 different pairs of nodes to be consi¬ 
dered); (c) obtaining the approximate value of the covariance between mean 
anomalies as 
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(4.16) 


r 

j 

k—1 k-1 k-1 k—1 rJ 

cov(igij,4g ki ) = £ E £ £ cov(4g"4gj tl ) 

■ s 0 isO rat 0 •= 0 

where m, n, r, and s are indexes that identify the elements of each pair of nodes 
in the subdivisions of a u and a kl . The covariances between point anomalies such 
as Ag“ and Ag kJ were obtained using a two term model for the degree variances 
of the anomalies: 

of** = 18. 3906—0.9943667 a+3 + 658.6132- j~ L -- 0.904894^ +2 

“ (n+100) (n+20)(n-2) 

This is the "2L" model of C. Jekeli (1978), and has the advantage that the value of 
the point covariance function can be computed using finite recursion formulas 
(Moritz, 1977). The same degree variances were used in (4.14) and (4.15). The 
number of point covariances being very large (k 4 ), a table with entries spaced 
at A = 0.05° intervals was created first, the required values being obtained by 
linear interpolation from this table. Numerical tests showed that k - 10 was 
large enough for both 5° and 1° blocks, because doubling this number resulted in 
a change of less than 0.2% in the values given by (4.16). Reducing the interval 
A^i from 0.05° to 0.005° had a negligible effect also, therefore the values obtained 
with (4.16) are probably accurate enough to test those given by (4.14). The only 
exception was the "polar” row, where the equal angular blocks are, in fact, tri¬ 
angles with a common vertex at the pole. Both with 5° and 1° blocks the dis¬ 
crepancies between (4.14) and (4.16) were large (more than 30%), regardless 
of how large a k , how small a At(i, or how big a N m9Y were chosen. The 
probable explanation is that the pole is given undue weight in (4.16) because it is 
treated as a whole row, instead of as a single point. For this reason, the 
numerical examples presented here stop at the row immediately below the pole. 

The first two tables show the covariances between 5° mean anomalies in the 
row between latitudes 0° and 5° (just above the equator) and in the tow between 
latitudes 80° and 85° (one below the pole). The error is at most 8%, though much 
less in most cases, and Njnax = 180 in each table. Under "Pellinen" one finds 
the values obtained using (4.15) with due regard for the change in block areas with 
latitude. While there is very good agreement near the equator, there is no resem¬ 
blance at all close to the pole to the other values listed. 


I 
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Table 4.1. Comparison between covariances of 5° mean anomalies Ag 
computed with expressions (4.14), (4.16) and (4.15), 
respectively. N max = 180. Row between 0° and 5° . 


Expression (4.14) 

Numerical 

Pellinen 

Block No. 

250.53 

253.77 

251.78 

1 

148.80 

149.46 

151.03 

2 

193.93 

93.95 

93.90 

3 

57.16 

57.15 

57.10 

4 

31.80 

31.76 

31.67 

5 

13.94 

13.93 

13.86 

6 

-18.09 

-18.09 

-18.08 

12 

9.12 

9.12 

9.11 

24 

-13.68 

-13.68 

-13.65 

36 


Table 4.2. Comparison between covariances of 5° mean anomalies Ag 
computed with expressions (4.14), (4.16) and (4.15), 
respectively. N max = 180. Row between 80° and 85° . 


Expression (4.14) 

Numerical 

Pellinen 

Block No. 

418.60 

437.23 

835.53 

1 

329.03 

318.27 

800.09 

2 

220.91 

229.17 

709.03 

3 

191.75 

196.40 

592.19 

4 

179.79 

179.70 


5 

165.58 

168.46 


6 

120.60 

123.73 


12 

69.89 

73.27 


24 

54.66 

58.01 

149.40 

36 


The next three tables show results for 1° blocks. Here the numerical 
method was conducted with the same A0 and k as in the case of the 5° grid. 
Results for rows between latitudes 0° and 1°, 45° and 46°, and 88° and 89° 
are shown. In tables 4.3 and 4.4 the discrepancies between (4.14) and (4.16) 
stay below 1%; this increases to about 5% near the pole (table 4.5). N ffaT is 
300 for the equatorial row, and rises to 400 from 45° on. As with 5° blocks, 
the "Pellinen" values are quite close to these of (4.14) and (4.16) near the 
equator, but become very different near the pole. 







Table 4.3. Comparison between covariances of 1° mean gravity 

anomalies computed with (4.14) and (4.16). 300. 

Row between 0° and 1°. 


Expression (4.14) 

Numerical 

Block No. 

849.68 

855.16 

1 

411.32 

410.30 

2 

220.35 

219.63 

3 

181.84 

181.50 

4 

163.33 

163.21 

5 

149.11 

149.15 

6 

1.07 

1.08 

31 

-16.79 

-16.79 

61 

1.99 

1.99 

91 

8.54 

8.53 

121 

-4.79 

-4.79 

151 

-14.39 

-14. 38 

181 


Table 4.4. Comparison between covariances of 1° mean gravity 

anomalies computed with (4.14) and (4.16). N max «* 400. 
Row between 45° and 46°. 


Expression (4.14) 

Numerical 

Block No. 

952.25 

959.17 

1 

531.15 

531.97 

2 

275.55 

274.65 

3 

210.81 

210.46 

4 

185.36 

185.57 

5 

170.82 

171.08 

6 

27.79 

27.79 

31 

-14.37 

-14.37 

61 

-18.01 

-17.01 

91 

-8.27 

-8.28 

121 

-1.01 

-2.92 

151 

1.39 

1.39 

181 
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Table 4. 5. Comparison between covariances of 1° mean gravity 
anomalies computed with (4.14), (4.16), and (4.15), 
respectively. N max = 400. Row between 88° and 89° . 


Expression (4.14) 

Numerical 

Pellinen 

Block No. 

1137.97 

1187.94 

1795.88 

1 

1135.79 

1184.15 

1795.02 

2 

1129.28 

1177.99 

1792.45 

3 

1118.55 

1155.08 

1788.18 

4 

1103.76 

1131.34 

1782.25 

5 

1085.15 

1102.85 

1774.68 

6 

408.59 

443.12 


31 

248.27 

257.88 


61 

207.19 

210.64 


91 

189.65 

192.55 


121 

181.93 

184.53 


151 

179.65 

179.86 

368.60 

181 


Calculations were carried out in the Ohio State University's AMD HAL 
470 V/6-n computer, using the FORTRAN H EXTENDED compiler, and double 
precision. The computing times for obtaining all 37 different covariances 
between the elements of a row of 5° mean anomalies, and all 181 covariances 
in a row of 1° mean anomalies, using (4.14) and (4.16) are listed in table 4. 6 
for comparison. The integrals of the Legendre functions required by (4.14), 
and the table of point anomalies covariances needed in (4.16), were precomputed 
and stored in core memory arrays, so the times given here do not include the 
determination of those auxiliary values. In most ordinary applications of these 
formulas, those quantities can be read from disk or tape whenever needed, 
because they are the same for a whole variety of problems. Clearly, using 
( 4 .14) can be orders of magnitude more efficient than using(4.16), while the 
accuracy is much the same. In fact, accuracy is probably better at the polar 
rows with (4.14) than with (4.16), because the latter seems to have problems 
handling triangular blocksi Finally, not all of the time-saving properties of 
(4.14) were exploited in the computer program used to calculate it, so there is 
3cope for some improvement in efficiency beyond that shown in the table. Notice 
the time saved in the 88° -89° row thanks to the neglect of terms in (4.14) that 
become too small near the poles, as explained in paragraph (4.2). 
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Table 4.6. Comparative efficiencies of algorithms based on expressions 
(4.14) and (4.16). 


Block Size 

Row 

N 

1N max 

Expression (4.14) 

Numerical Quadratures 

5° 

u* 

1 

o 

o 

— 

0.1 sec. 

23.3 sec. 

5° 

80°-85° 


0.1 sec. 

23.3 sec. 

1° 

0°- 1° 

w 

0.35 sec. 

113.3 sec. 

1° 

45°-46° 


0.43 sec. 

113.3 sec. 

1° 

88°-89° 

400 

0.06 sec. 

113.3 sec. 


4.4. Covariances Between Mean Values and Point Values 


The prediction by least squares collocation of mean values from point 
values, or vice versa, requires finding the corresponding ''mixed” covariances. 
In such case, formula (4.2) becomes 


cov(u n , v(0,A)) = ~~ f cov(u(P'), v(P)) da 

J<T u 


(4.17) 


where P' = (0', X') e <T U , P = (0,X), v(P) is the "point” value of v, and u u 
is the average over the i, j block of the grid, as before. Following a similar 
reasoning, one arrives at a formula that corresponds to (4.14), except that only 
one area integral has to be considered. The new expression is: 


X N,ax 


u y X 

COV(U U ,V(0,X)) - ~ „£ 0 [ h £o b 5.( 2n +1) ■RijSNb +d( coS ®) ) 


N nax 


1 

— _ r 
+ E L b ( 2 °" 1 ) ( I n^Nh-a,l P n ,3Nh- a (cos0)J [(A(m) cos m(Xj - X) 

+ B(m) sinm(Xj-X)J (4.18) 


where 


A(m) 


= I 


AX 


if m = 0 


((cosmAX-l)/m otherwise 


B(m) 


I (sim: 


0 if m = 0 
(sinmAX)/m otherwise 


and the ( I n>n>1 , the c“ ,v , and the N, K, and N max areas in (4.14). Expression 
(4.18) assumes that the area means belong to a grid with constant AX. If the point 
values are also on a grid, and if this grid is congruent with that of the mean anom¬ 
alies, implementation of (4. 18) is quite efficient. In fact, the speed of a good 
algorithm for doing this should be much the same as that of one for implementing 
(4.14). On the other hand, computing the same covariances by numerical quadra¬ 
tures is k s times faster in this "mixed" case than it was in the previous one, be¬ 
cause there is only one area integration involved. Assuming k = 10, as in the 
previous examples, then expression (4.18) should be only 2-3 times faster than 
numerical quadratures for 5° mean values, and from 4 to 15 times faster in the 
case of 1° averages. 
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1 


Both the normalized Legendre function and their integrals are needed for 
(4.18) . They can be precomputed and stored on disk or tape until needed, in 
this study the following recursive relationships were used to generate their 
values: 

p— v - l( cose, - ) ]*c M eiu,.„ ( co s e> 

r(2n - l)(n + m - 3)(n - m - 

“L (2n-;5)(n + m-2)(n-m) J p n-3,»-i( cos 0 )» m^n (4.19a) 
p.-i,.-i(oose> = sine (cos9) 

J^ a p.-i,.-i(coa0)sIn9(J8 - '-J[-ff^ 2 + ° m ".\]'sin s ep._ a> ._,(cose) 


(4.19b) 


, (n-3) r (2n-l)(D + m-3)(n-m-ln ^r 6 a_ « . „ ^ 

n L (2n-5)(n + m-2)(n-m) J J e p »-3,«>-i( cos 0) sin 0 d0 (4.20a) 


Ga 

6i 


f e s- - _ _ 1 r (2n-1) „ 

J 01 p n-i,„-i(cose)sin6 d0 = 2^|_(n-l)(n-2)J sin 9 p »- 3 ,n- 3 (cos 0) 


0a 


n 

2nL 


- l)(2n - 1)(2 n- 3H $ p9a 


(n-2) 


sin0 d0 


(4.20b) 


if 0 is not very small, or 
• 0 
0! 


i 02 


J eap n-i,„<.i(cos0) sin 0 d0 = ~ • • 3 1 ^ ; »+i 0 

J 0j. l_(2n-2)(2n-4)...4J sm 

r 1 + 1 sin a 8 1*3 sin 4 6 1-3-5 sin^ -| 

Ln + 1 2 n + 3 + 2-4 n + 5 + 2-4-6 n + 7 + "’J 

if 0 is very small. Here y( 0) |!® = y( 0 a ) - y( 0 X ). The recursive form 
the integrals of the Legendre functions were derived by M.K. Paul (1978); 


(4.20c) 


formulas for 

^ _ _ __ __ < v tlic 

author has been fortunate enough to have available a FORTRAN subroutine pro¬ 
grammed by Paul, and kindly sent by him to Professor R.H. Rapp of the Depart¬ 
ment of Geodetic Science at O.S. U. The results reported here have been made 
possible by, and bear witness to, the great numerical stability of Paul's for¬ 
mulas. 
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5. Conclusions 


The relationships between spherical harmonic series and Fourier series, 
coupled to the symmetries of spherical grids, permit the development of ef¬ 
ficient algorithms for numerical analysis of data regularly sampled on the sphere. 

The algorithms presented In section 1 for implementing numerical quad¬ 
ratures are efficient enough to allow the analysis of 64800 1° x 1° mean values 
through degree 180 in less then 20 seconds, and the summation of the 90000 
terms of a harmonic series complete through degree 300, on a full 1° x 1° 
equal angular grid, in less than 1 minute. The analysis of large global data 
sets to a very high resolution is a relatively trivial operation with modem 
digital computers. 

The principles and ideas behind the optimal estimators of section 2 pro¬ 
vide a rational basis for the study of linear techniques for spherical harmonic 
analysis, both optimal and non-optimal. The error measure introduced in this 
section is shown to be a very reasonable way of evaluating the estimation error, 
as illustrated by the results listed in table (3.1). 

The optimal estimators themselves are reasonably easy to compute and 
use, particularly when they are of the quadratures type, which happens under 
the fairly general conditions discussed in section 2. Even when such conditions 
are not present, as in the case of scattered data, the problem still has a 
structure strong enough to allow efficient algorithms for creating and inverting 
the normal matrix. 

The separation of the problem of estimating the coefficients (by least 
squares collocation or by least squares adjustment) according to the order m of 
the coefficients, allows both for efficiency and for numerical stability. Even 
if the total number of unknown coefficients is very large, the largest matrix 
to be inverted is of dimension N , instead of 0(N 2 ), as it would be if the prob¬ 
lem could not be separated In this way. 

All the algorithms presented here, when the grid is complete and regular, 
are well suited to parallel processing. 

In the case of full grids of mean values with even noise, the results of 
section 3 suggest that the optimal "collocation" estimator can be approximated 
very closely by a much simpler quadratures-type formula, the "composite for¬ 
mula" (3.9). The search for simple, near-optimal estimators is just as impor¬ 
tant, from a practical point of view, as the search for efficient algorithms for 
obtaining the optimal estimators themselves. This is a topic that certainly 
deserves Ihrther research. 

The methods for creating and Inverting the normal matrix, that make 
possible to find optimal estimators for large data sets, have application outside 
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spherical harmonic analysis, in all areas of estimation, filtering and prediction 
on the sphere. This has been the subject of a previous report (Colombo, 1979). 

It must be added that the principles presented here can be generalized to bodies 
of revolution other than the sphere; to the case where the data are not homo¬ 
geneous (i.e., a mixture of, say, gravity anomalies, satellite altimetry, etc.); 
to the case where the coefficients to be estimated are not those of the signal as 
given, but of some more or less complex linear transformation of this signal, 
satellite-satellite tracking data being a good example. In fact, the author is 
at present considering the error analysis of the determination, by least squares 
collocation, of the potential coefficients from satellite-satellite tracking, fol¬ 
lowing some of the principles of section 2. This will be the subject of another 
report. 

The method developed in section 4 to calculate the covariance between 
two area means without employing cumbersome numerical integrations is of 
interest, not only in spherical harmonic analysis, but more generally in fil¬ 
tering, prediction and estimation from mean values on the sphere. 

The computer programs described and listed in Appendix B should help 
the interested reader to implement some of the techniques discussed in this 
report. The author sincerely hopes that this will be done by workers concerned 
with improving and further developing such methods. 

Above all, the author hopes that he has conveyed to those who had read 
this far, the idea that the detailed analysis of vezy large sets of global, reg¬ 
ularly sampled data can be done within the computing resources available today 
to most scientists who work at universities and research institutions every¬ 
where. The processing, be it by numerical quadratures, or by simultaneous 
adjustment, of "all the data in the world" is not a fanciful thought, but a practical 
possibility. 
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Appendix A. Term by Term Integration of Formula (4.5), and Rearrangement 

of Terms in Formula (4.9) to Arrive at Formulas (4.11-a) and (4.14) 


The area integrals in expression (4.5) can be split into integrals in 0 and 
X, and the summation theorem for fully normalized spherical harmonics can be 
used to replace P a (^ PP ') with an equivalent expression in the general term of the 
series; 


_ _ i (•9i + Afl 

cov(u 1J ,v !tl ) * -—— | sin0d0, 

A *i A w J 9 t J 9t! 


I* r Xj + AX » Xi + AX ® 

" A* " dX i dX’ £ c'i’ v . 


% 


£ P aa (cos 9 ) Pa, (cos 0' ) cos m( X - X') (2 n + If 1 


(A.l) 


B =5 0 


Starting with (A.l), the proof proceeds in four steps, each justifying in turn 
the taking of one of the four integrals inside the summation symbols. Each time, 
three theorems are invoked: 


The first is the "M-Test" theorem, due to Welertrass (see, for instance, Carslaw, 
(1950)): A 

The series f(x) = o 4 i 0 u (x) n will converge uniformly in a^x^b 
if there is a convergent series of positive constants M 0 , M T ,..., M„,... 
such that, for all x in a^x*b, |u(x) n |£M n for every positive 
integer n. 

The second theorem is (also according to Carslaw): 

00 

If the general term u(x)„ of the series JC 0 u(x) n is continuous, 
and if the series converges uniformly to some function f (x) in 
the Interval asx^b, then 

p b p b « rnPb 

;■ f(x)dx = £u(x) a dx = £ u(x) a dx 

This is a sufficient condition for term by term integration. The third theorem 
is the mean value theorem for integrals 

If f(x) is analytic in a«xsb then J^f(x)dx = (b-a)f(c) 
for some c such that a^csb, 


Proof : The series in the left hand side of (A.l), if all variables but X' are kept 
fixed, is uniformly convergent in the interval Xi £ X' 5 X: + AX because 

n 

(2n + if'cn ,r E 0 P„ (cos 9) P n , ( COS 0') cos m( X - X’) = c,“’ T P a (') 

and |c a ,Y P a (&p>) j s c a ,r because max | P,(i p )! * 1 (the argument of P a is real) for 

0 s * £ tt 
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ail n, while £ 0 c B * v * cov(u(P),v(P)) ia always finite (equals the value 
of the ’^>oint" covariance when = 0 ). The "M-Test" condition is satisfied 
so the series converges uniformly; therefore, term by term integration with 
respect to X’ is valid: 


cov(u 1Jt v kl ) » — j sin0d0j sln0* d0'j dX • 

t« hi 9j 0 k Xj 


©W-A0 


r Xj+AX 


• T c»“ ,y £ £ a ,(coe§) ^.(cosS') AX cos m(X-X') dX(2n + l) -1 
n »0 a> 0 J X 1 


Applying the mean value theorem to the last express bn: 


■j « 0 i + A 9 (> 01(00 

cov(Un, v kl ) = -—— sinSdS sinfl'd0' dX • 

J A n A k i J 0t Je k 

oo a 

• £ Ca’ r £ AX Pa«(cos0) P B ,(cos0') cos m(X-Xq) (2n + l) -1 

B=0 *»0 

where Xj ^ Aq ^ Xi + AX . Removing the common factor AX from the summation: 


e k + Ae 


,Xj + AX 


*\ «9i + A 9 <ok-«0 »A. 

cov(U!,v kl ) = — — I sin0d0 | 8in0' d0* 

J e t J e k J Xi 


0 k + A9 


,Xj + AX 


■ £ c»’ T £ Pa.(cose) 2 a .(cos0’) cos m(X- x,)(2n + l)" 1 (A.2) 

a = 0 r < = 0 


The general term in the partially integrated series is 

n 

(2n + l) 1 c a u ’ T £ P tt ,(cos0) P M (COS0') cos m(X-X)) = c„ u » T P n ((M 

■ =sO 

where Q = (0',Xq); now |ci ,T P B ( $?q) I s c B ’ r for all n, and for all X in the 
interval Xj ^ X s Xj + AX , so the "M-Test" is satisfied again and the series is 
uniformly convergent, and thus integrable, with respect to X . Therefore 

»\3 _9j + A9 9ic + A9 <* 

cov(u t4 ,v kl ) = —— J sin0d0 J sin0'd0 J] c u / • 

a 

• £ Pam(COS0) P B1 (COS0') cos m(X* - Xq) (2n+l) _1 (A. 3) 

• = 0 

where Xj £ X* * Xj+AX. Once more the general term of the twice integrated 
series satisfies the "M-Test", because 

v n 

|(2n + lf 1 Ca’ Yj Pb»(cos0) P m (cos0’) cos m(X«-Xq) | * |c B ,r p B (^)| * c B ,r 

•*o 

(where R s (0,X*)) for all n and, in particular, for 0* in the interval 
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0 k £ 0' s 0 k + A0. Inconsequence 


cov(u t <, v kl ) = 


4A a sin$j <40 Z *® 1 a Jfi O “*■' 
——— stoSdS £_c. 


u kl 
c 


“e* 


c ■ 0 


• £ Pn,(cos0) P„(CO80 S ) COS m( A*- A„) (2n + l) _i (A.4) 


•=o 


Finally, 


|(2n + l) -1 c n u,Y £ P„(cos0) P M (cos0s) cosm(A»-A Q )| £ |c n v FU&*) | £ c a ’ v 
*»0 

(where X = (©s,Aq)) for all n and, in particular, for 0 in the interval 
Qi £ 0 £ 0j + A 8 so the last integral can be put inside of the summations, and the 
proof of the term by term integration of (4.5) is complete: 


cov(u u ,v kl ) -- h£ Cn* v £ ^.(cosOr) P„(co8 0 s ) cosm(A(*-A g ) (A.5) 

B S3 0 *~0 

where 0, £ 0r £ e t + A0 and h = AA a sin 03 £0 sin 6, AO/A^.A^ 

The general term in (A. 5) satisfies the "M-Test": 

a 

u v - u V 

|z n |-|(2n + l)- l c a ’ £ P».( cos0r) 7a»(cos0 s ) cos m(A»-Aq) I S c B ’ 

■ sO 


Since the series a £ 0 Cg ,T converges, any series of positive terms |z„| satisfying 
| z„ | < c B must converge also: the ' t M-Test" condition implies the absolute con¬ 
vergence of a g o z n . Absolutely convergent series can have the order of their 
terms changed arbitrarily, without changing the value of their limit sums. This 
justifies the reordering of (4,9) that leads to expression (4.11-a). 
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Appendix B: Computer Programs (Descriptions and Listings) 


This Appendix contains the description and listing of each of the major 
subroutines developed in the course of this research, together with their own 
auxiliary routines. The listings of the main subroutines contain explanatory 
comments that the author hopes will be of help to those who may use them. 

The description accompanying each main program defines arguments, gives 
the dimensions of the memory arrays required, mentions the relevant for¬ 
mulas from preceeding sections, and gives a brief explanation of the various 
segments of the program. 

B. 1 General Programming Considerations 

Only one-dimensional arrays are used in the software described here. 

The language used for all of them Is FORTRAN IV ; some subroutines from 
the International Mathematical and Statistical Libraries Inc. (IMSL) are called 
by the main subroutines. Implicit DO loops in READ or WRITE statements have 
been avoided as much as possible, because their execution may be rather slow, 
depending on the compiler; instead, subroutines FREAD, FWRITE and REWIND 
are used for all imput/output operations involving large files on tape or disk, 
in some of the subroutines. All operations involving real arithmetic have been co 
ded in double precision ( 8 bytes, or 32 bits), which is equivalent to retaining 
the first 7 significant figures in all arithmetic operations. 


The arrays containing the associated Legendre functions or their integrals, 
as the case may be, are arranged first bv degree, and then by order : 00, 10, 
11, 20, 21, 22, . . . (N aaa , N„, x ).The ££,are arranged accordingly, always in 
two separated arrays: one for the c£n another for the C n \ . In order to 
get the value of the element "nnt' from one of this arrays, the following formula 
is used; 


k = £n(n + l)+m+l 

where k is the position of this element in the one-dimensional array. When 
the elements are recovered sequentially from the beginning (00), the following 
type of DO loop is used; 


KOUNT «1 

Do XX N1 = 1, N aax 

DO XX Ml = Nl, N aax 

LEGEND (KOUNT) » ARRAY (KOUNT) **2 

XX KOUNT = KOUNT + 1 

where, in this particular example, the nm (n=Nl-l, m=Ml-l) element in array L 
LEGEND is equated to the square oi the nm element in ARRAY. Avoidance of two- 
dimensional memory arrays results in considerable improvements in efficiency. 
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B.2 Subroutine SSYNTH 


This subroutine computes the sum of a spherical 
harmonic series complete to degree and order NMAX at each one of the 2N 2 
points or blocks in an equal angular grid. The subroutine can calculate point 
values (IFLAG = 0) or area means (IFLAG = 1) . The number of rows or paral¬ 
lels N (Nyquist frequency) must be even . Subroutines FFTP from the IMSL Double 
Precision Library is used to calculate the sum of the series along rows by means 
of the Mix Radix Fast Fourier Transform algorithm. 

The procedure used is that described in paragraphs (1.6) and (1.7). The 
symmetry of the grid with respect to the equator, and the corresponding even- 
odd symmetry of the values of the Legendre functions or their integrals, (I.e., 
the x“* of section 1) are exploited. The values of those functions, or of their 
integrals, are read from mass storage (disk or tape) into array ROW , in the 
order described in the previous paragraph. All the values for one latitude, or 
tow;*are read at once, so the dimension of ROW is £(N i>x + 1)(N,. X + 2). All the coef¬ 
ficients Cf, are also stored in core, the cor esponding RCNM and RSNM ar¬ 
rays (for C‘„ and 5*, respectively) have the same dimension (the are 

included, though they are all zero). The output consists of 2N a values in ar¬ 
ray DATA . This array is organized in rows, from North Pole to South Pole. 

The rows, of 2N points or blocks each, have their values written consecutively. 

The following is the list of arrays, and their dimensions: 


NAME 

DIMENSION 

TYPE 

ROW 

RD4(N„ x * 1)(N,. X + 2) 

REAL *8 

RCNM 

RD 

Tf 

RSNM 

RD 

Tf 

X 

RD 

INTEGER * 2 

DATA 

2N a (N = 180/BLOCK) 

REAL *8 

CR1 

N + 1 

Tf 

CR2 

N + 1 

Tf 

SRI 

N + 1 

»f 

SR 2 

N + 1 

Tf 

AM 

N + 1 

tf 

BM 

N + 1 

ft 

F AUX1 

4N 

Tf 

F AUX2 

4N 

•T 

F IWK 

see IMSL Handbook 

INTEGER * 4 

F LL 

IT ?T Tf 

LOGICAL *4 

F A 

ft ft r* 

REAL *8 

IV 

N + 1 

tf 


("F" designates those arrays required by the IMSL subroutine FFTP). to ad¬ 
dition, the sice (in degrees) of the blocks Is defined by BLOCK; WPP * i(N,» x + 1) 
(K m ,+ 2); IU is the number of the unit (disk, tape) from which the Legendre 
functions or their Integrals are to be read. 

Array X contains information on whether a given x?" is even or odd; 
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arrays CR1 , SRI ,and CR2 , SR2 , respectively, contain the Fourier coef¬ 
ficients a, (ai) and of rows i and N-i-1 ; arrays AM and BM 

contain the values of A(m) and B(m) , as defined in expression (1. 7); the 
auxiliary array IV contains the numbers £n(n + 1) needed to locate individual 
elements within arrays ROW , RCNM , and RSNM , when they are not addressed 
sequentially. 

The comments in the listing are probably enough to understand most of it 
on close inspection; one point however may be worth explaining further: the 
"aliasing" of the Fourier coefficients has been incorporated to take care of the 
case when N„„ > N - 1 . In such situation the a, (ai) and S^b 1 ,) become 
aliased, as Fourier coefficients must, and it is their aliased values that the FFT 
subroutine requires to compute the values of the spherical harmonic expansion 
along parallels. The formulas for the aliased coefficients are 

KM 

a 1 . 3 a, + 3C (°d«N + «*!*-») (Bl-a) 

- si + r (S***n - An-i) (b 2-b) 

where KM is a large enough integer. A similar expression applies to a, 
and to b, . 

The arrangement of the output in latitude corresponds, in the case of area means, 
to the intervals on which the P B , are integrated; for point values, It is de¬ 
fined by the latitudes 9, at which the P B , (cos 9 t ) have been precomputed. As 
regards longitude, the grid starts from the zero meridian used for defining 
the coefficients. In the case of point values it Is usual to compute all values at 
the center of each block. To do this, the F„ . must be precomputed at the lati¬ 
tudes of the center points, while the longitudes are taken care by modifying the 
coefficients as follows 

(?»’, = C n , cosm y- + S" al sinm (B 2-a) 

S a '« = S’q, cos m -C al sinm^ (B 2-b) 

This is equivalent to rotating the grid eastwards from the zero meridian by ^ . 


B.3 Subroutine HARMIN 

This subroutine implements either the algorithm of paragraph (1.5) 
for the harmonic analysis of area means, or that of paragraph (1.7) for the 
analysis of point values. 

The subroutine calls IMSL's FFCSIN to calculate the a 1 , , b[ or the 
by means of the Fast Fourier Transform (Mix Radix) algorithm, it 
also calls subroutine QUADFS , that returns in array A the de-smoothing 
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factors . QUADFS calls LEGPOL , a subroutine that computes the Legendre 
polynomials up to degree NN + 1 needed for the fT^ in (4 . 


The data is arranged as in SSYNT1I , in array DATA , before the subroutine 
is called. Afterwaxds, the contents of DATA are destroyed, as the a a , b B 
or the a i , 0 B are formed in place of them, row by row, by FFCS1N . The 
resulting coefficients' estimates are put into arrays RCNM and RSNM , in the 
same order as for HARMIN . The other arrays, with the exception of A , are 
as in SSYNTH . The same is true of the scalars, with the exception of NN . 

NN is the highest degree and order to be estimated. NDD is the total number 
of Legendre functions,or their integrals, lobe read from unit IU , per . 

This number is (NN+1)+(NN+2)/2. A is a REAL *8 array of dimension NN + 1 . 
The dimension in QUADFS allows for a maximum NN = 300; for larger solutions, 
the dimensions there and in LEGPOL must be increased accordingly. 

In the case of point data, the estimated coefficients are computed using a 
center point formula that assumes that the data are situated at the centers of the 
blocks; the resulting coefficients are referred, nonetheless, to a grid starting 
at the zero meridian (the "rotation" of the coefficients takes place between 
statements 0071 and 0073, when IFLAG = 0) . When IFLAG=1 , the area means 
formula (1.30) is computed; the xV = MnJj Ba (8) sin 9 dp , and the /i B are 
those produced by QUADFS, as already mentioned. The integrals of the Legendre 
functions are read from unit "IU", as in SSYNTH (same format), and the size 
of the blocks is specified by BLOCK (in degrees). The version of QUADFS 
listed here implements the "composite" estimator of paragraph (3.3). If 
another Is desired, this can be achieved simply by replacing lines 0021 through 
0024 in QUADFS. 


B.4 Subroutine NORMAL 

oc 

This subroutine creates the optimal estimatox*s for the C aB based on the 
formation and inversion of the R(m) matrices described in paragraph (2.10) 

The algorithm exploits the fact that (C zi + D) is a block matrix of Toepiitz 
circulant sub-matrices. This subroutine is meant only for mean values. 

The grid is as in SSYNTH and HARMIN. The symmetry with respect to the 
the equator is only partly exploited* matrix D may not be persym- 
metric, so the total matrix (C„ + D) may not be so either. C zz however, is 
always persymmetric, and this is taken into account to save computing and stor¬ 
age. A general diagonal matrix D corresponds to a rather broad class of 
actual problems,such as the analysis of the 5° x 5° real gravity anomalies 
described in paragraph (3.4). 

This subroutine requires four input/output units: 8 (read only) contains 
the values of the integrals of the Legendre functions, row by row, arranged as in 
SSYNTH or HARMIN; 10 contains the right hand sides of the "reduced normals" 
.lira = R( m ) 2£ n * (expression (2.58)); 15 contains the R(m) matrices, ordered by 
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increasing m , stored in vector array form column by column; 30 stores the 
Xi* of the optimal quadratu res -type estimator. The xi* are stored from N. 
Pole to S. Pole, and according to nm , as the Legendre integrals and the coef¬ 
ficients. In some circumstances the grid may be geocentric rather than geodetic 
and a change of coordinates might be desirable: this can be achieved by setting 
the parameter IGEO to 2 . The flattening assumed for this transformation is 
F = 1/298.257. 


After the R(m) matrices have been created, they are inverted by IMSL 
subroutine LUDECP , that performs a Choleskii factorization. IMSL sub¬ 
routine LUELMP solves the equations resulting in the Xi" ; if during the in¬ 
version LUDECP detects an ill-conditioned (or a singular) matrix, the solution 
part is avoided, and a set of null Xi* is stored for that particular m . As 
an additional check for the stability of the solution, the relative residuals. 

N-l 

r v » t 

r ' “ *-° where v =[v 0 , . . . v^J is 



v = k*" - R(m) x aB (computed) 


(B 3) 


are computed and printed. In all the cases studied here these residuals indicated an 
agreement of at least 9 significant figures. To improve the stability of the 
solution, a regularizing constant REGUL is added to the diagonal elements of 
the R(m) (Paragraph (3.3)). 


Arrays PN , SS , and FC contain the propagated noise, sampling, and 
total error measure (variance) per degree. W contains the averaged row 
variances (expression (2.43)) arranged from North to South. 


The scalar arguments, NMAX , NN , DGRID , IGEO , REGUL , NRUN , and 
NC 2 , are described in the comments inserted between statements 0006 and 0010. 
The arrays are as follows: 


NAME DIMENSION 


TYPE 


ROWP ■§(NMAX +1)(NMAX+2) REAL *8 

ROWQ " " 

RHS i-NN (NN + 1) N (N = Nyquist freq.) " 

S (NN + 1) N " 

A kl(N - 1)+N " 

UL " " 

W N " 

DVAR NMAX +1 " 

FC NN + 1 " 

PN " " 

SS 


Arrays ANMPQ , FF, XO , B , X , BT , FINMP ( all REAL* 8), and IDD (REAL* 4), 
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all dimensioned 200 or 400 in the subroutine itself, are large enough for 
problems where N < 200 . For finer grids, the size of these arrays should be 
increased in die same proportion as that of N . 

The R(m) matrices are formed according to expression (2.63) . Subroutine 
FUR computes the "aliased" Fourier coefficients of the covariance functions 
that are, in fact, the elements of the R(m) , scaled by N or 2N , depending on 
m . Common MM and array MT are part of a logic set up to ensure that the 
Fourier coefficients are not computed more than once each. 

Subroutine ANALYS uses the Xi* stored in unit 30 to analyze the data 
in array DATA . The rc, 1 , are formed in place, as in HARMIN , so the 
original values in DATA are destroyed. Arrays CR , SR , CAA , CBB , SAA , 
and SBB have all the same description as CR1 , SRI , etc., in HARMIN . 

IMSL subroutine FFCSIN (double precision) is used to obtain the a, 1 , /3, 1 . The 
reason why ANALYS is used instead of HARMIN , is the arrangement of the 
X*" "columnwise", or by increasing latitudes, rather than "rowise" (i.e., all 
the Xi" for the same i stored together) as HARMIN would require. 

The listings of FUR , ANALYS , and those of the fast input/output sub¬ 
routines FREAD , FWRITE , REWIND , are given after that of NORMAL . 

The input/output subroutines have dummy arguments, because originally NORMAL 
was written to work with certain subroutines available at O.S. U. that may not 
be in the software libraries of other institutions. 


B. 5 Subroutine NORMAX 


A modified version of NORMAL , this subroutine was created to compute 
the variance of the estimation errors in ordinary quadratures formulas according 
to the theory in section 2 . Essentially, it computes 


'€ » 
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°n~2 Y. ^ £««fl 
,=o a-o 


» i- ft T 

+ I T (f“.) (C„ + D) 
■ =o Of=o 



by forming and using the R(m) . Since no inversion or solution of the normal 
equations is required, the corresponding segments have been removed from 
NORMAL , and a new final segment added for the computation of the various 
accuracies. 


The theory behind the calculation of the o^ a using the R(m) matrices Is 
as follows*. 


In the case of ordinary formulas of the type (see expression (1.7)) 


- 124 - 













IT Cl RELEASE 2.1 KMUUX DATE • 8*243 12/29/84 PACE M8S 



UUUU UUU WWW wwwww 



|t:8S|HI 



127 









IV Cl ■«■«« >.• MHUIAX DATE • 8*342 12/24/34 PACE M44 



128 

























6“ = (f n “) T m 

N_1 3N-1 /. _ 

= 1/ Y »“M dff mo 

1-0 

= r fxrt cosmjAX + | sinmjAxlmjj 

i—o L j=o*dB(m)J *• A(m)-* J 

the estimator vector is a combination of a "sine" and a "cosine" vector of 
"frequency" m (terminology introduced in paragraph (2. 10)). The product of 
such vector by (C zz + D) is, because of the structure of this matrix, another 
combination of a "sine" and a "cosine" vector of the same "frequency". From 
the pi’operties of the normal matrix follows that 

cos m j A X + 

B(m) \ 

. W sinm ) 4X • 

where, calling 

y n °-[ Vg\ i/“; . . . z/ N n ~\] T = R(m) [xo>» XT> • • • x“i) T 


and 

F(m) 

( 2N AX 2 if m = 0 or m = N 



14N (1- cos m AX) if m ^ 0 , m ^ N 


is 

& 

n-i 



(C„ +D) f® = F(m) £ XT*?' 

(B4-a) 


(C„+D)J“ = 


2N if m = 0 
or m = N 
N if ra / 0 


jA(m)\ 
Wm) ' 


r a 

Regarding the scalar product 2c oa0!fI _f n , , it is easy to show that 

2 £*ot.z f“ - Z 2^TlX (X?#)S F(m) (D4 ' b) 

Expressions (B4-a) and (B4-b) are implemented by NORMAX to obtain the 
error variances. This subroutine also uses FUR . Subroutine QUADRS is 
also called, to obtain the de-smoothing factors. In the case listed here, this 
factor is = gbr . Array QUADS has been added to the list of arguments, 
and it contains the after the call to QUADFS . 


B. 6 Subroutine LEGFDN 

This subroutine computes the normalized Legendre functions and, if so 
desired, their derivatives at a given 0 t . All values correspond to the same 
order M ; if more than one order at a time is needed, a DO loop, where the 
subroutine is called once for each order, can be set up. The subroutine is 
based on formulas (4.19 a-b) and (1.38 a-b). The use of this subroutine is 
explained by the comments inserted In Uic listing. The stability of Die recursive 
formulas was tested by computing P n , (cosO) and (dP»,/dG)(cosO) for m - 350 
and 350 -■ n • 400 , 2.5° sqs 90°, at 5° intervals. Calculations were done 
first in double precision (8 bvtc words) and then repeated in extended precision 


- 131 - 







r 


i 

s 

a. 


s 

\ 

8 

£ 


(i 


S §? 


*1 

£f 

ii 

2b 


i 1= 

ci -A 


s 


. ?s Sslli 

e i§ =•«:; 
1 gS*53 

*S 

a«s:"i5Pp:p 

W§5-5«5 

88St£*t: 


SS 


s 



WWW 


> 


f !fi fill!! illiliifiiiilii ill 


133 





(16 byte words). Both sets of results agreed with each other to better than 
10 significant figures. 


B. 7 Subroutine NVAR 

This subroutine computes the degree variances of the gravity anomalies 
(point values),up to degree 100, according to the coefficients obtained by 
R. Rapp from a complete, equal angular set of mean 1° x 1° anomalies, as mentioned 
in par. (3.1). Above n = 100 , the subroutine uses a two-term model to calcu¬ 
late the 9 |(Ag) . The resulting degree variances are stored in array DVAR . 

The first element in DVAR is uf (Ag) . 
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B. 8 SUBROUTINE COVBLK 


Subroutine COVBLK calculates covariances between area means according to expres¬ 
sion (4.14). This subroutine is listed in the following pages. The arguments are 
explained in the comments at the beginning of the listing. In addition, the 
following things have to be born in mind: the dimension of the array DVAR 
is N max ; the dimension of both RINS and RINN is iKNj^ax- 1 - 2 )»* the 

dimension of COVS is 360/BLOCK. The values returned in DVAR are the 
original degree variances c“’ , each divided by 2n + l. If LB is less than 
360/BLOCK, the LB + 1, LB+ 2, ... ,(360/BLOCK - LB -1) elements in COVS 
are returned as zeroes, the remainder contain the first (and last) LB covari¬ 
ances. The dimension of F is 180/BLOCK (Nyquist frequency). To use this 
subroutine with Nmax > 400, the arrays FF and IDD (whose dimension should 
be no less than Nmax)i should be redimensioced. 

The subroutine does not take advantage of the ''aliasing" of Fourier coeffi¬ 
cients built into expression (4.14). Implementing this aspect should lead to 
some additional improvement in efficiency. The Fourier series is computed, 
orcn th~ .-.-efficients have been determined, by multiplying each coefficient by 
the uorrospooning cosine of mAj and adding the products together. The values 
of cos m/v f are computed using the following recursive formula: 

coa eiAj * 2 cosAj cos(m-l)Xj - cos(m-2Uj 

which avoids repeated calculation of the FORTRAN COS function (only cos A. is 
required to start the recursion). Actual calculation of the Fourier series requires 
about 0.04 seconds in the most time consuming case: the grid of 1° blocks. The 
greater part of the time taken by this subroutine goes into finding the Fourier 
coefficients of the mean value covariances. For this reason, there is not much 
difference between computing all covariances in a certain row, or just a few of 
them, using this procedure. 
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